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1.       INTRODUCTION 

In  a  recent  report  [1]  a  mixed  numerical/analytical  approach  to  the  computation  of  a  ring- 
symmetric  spacecraft  exhaust  plume  was  presented.  The  numerical  scheme  had  been  implemented  in 
a  code  named  "JET"  which  is  capable  of  generating  whole-plume  flow  fields,  while  the  analytic 
approximation  is  restricted  to  the  ring- symmetric  centered  rarefaction  waves  (CRW)  that  flank  the 
plume.  The  present  report  is  intended  to  serve  as  a  supplement  to  [1]  in  providing  details  on  the 
computational  scheme  and  the  code  JET. 

The  spacecraft  exhaust  flow  (Fig.  1  of  [1]  )  is  idealized  as  a  ring- symmetric  steady  isentropic 
expansion  of  an  ideal  gas.  The  nozzle  lips  are  assumed  sharp;  the  supersonic  flow  from  the  exit 
surface  of  the  ring-nozzle  is  assumed  uniform,  and  the  background  is  considered  to  be  perfect 
vacuum. 

The  standard  scheme  for  computing  such  idealized  ring-plumes  is  the  classical  (direct)  method  of 
characteristics  [2]  .  At  a  preliminary  phase  of  the  present  laser  exhaust  study,  a  code  AXS  YM  [3]  was 
written  for  computing  ring-plumes  using  this  method.  A  notorious  shortcoming  of  the  direct  method 
of  characteristics  is  that  the  solution  grid  is  highly  irregular,  being  formed  by  the  (oblique) 
intersection  of  the  C  and  the  C  ~  families  of  characteristic  lines.  We  first  encountered  a  difficulty 
with  this  grid  while  seeking  a  scheme  for  integrating  the  molecular  opacity  along  a  straight  line  [1]  . 
This  computation  would  have  required  rather  complex  coding  for  the  geometry  of  intersection 
between  a  straight  line  and  an  irregular  grid.  It  seemed  preferable  to  opt  for  a  computation  scheme 
that  would  produce  a  more  regular  grid,  even  at  the  expense  of  some  loss  of  accuracy.  Such  scheme 
is  the  inverse  marching  method  of  characteristics  [4]  . 

Generally  the  marching  in  this  type  of  scheme  is  in  the  downstream  direction,  i.e.,  the  y  direction 
in  our  case.  Grid  points  are  located  on  a  succession  of  constant-y  rows,  thereby  introducing  a 
measure  of  regularity  in  the  solution  grid.  Just  two  rows  have  to  be  stored  in  the  computer  core 
memory  -  the  "old"  row  and  the  "new"  row,  whereas  in  the  direct  method  o[  characteristics  whole 
grid-image  matrices  are  required  to  reside  simultaneously  in  core  memory. 

The  first  version  of  the  JET  code  was  based  on  the  inverse  marching  scheme  given  by  Zucrow  and 
Hoffman  (Section  12-5  in  [4]  ),  where  the  flow  variables  were  the  two  cartesian  velocity  components. 
The  computation  seemed  accurate  everywhere,  except  within  the  centered  rarefaction  wave  (CRW). 
In  an  attempt  to  replicate  a  planar  CRW  (Prandtl- Meyer  flow),  the  numerical  solution  exhibited  an 


instability  :  Mach  number  increased  along  the  (low  pressure)  boundary  characteristic  line,  rather  than 
remain  constant. 

A  qualitative  explanation  for  this  instability  is  the  following.  Flow  gradients  in  a  CRW  are 
inversely  proportional  to  distance  from  the  corner,  so  that  the  inverse  marching  scheme  gives  rise  to 
an  amplification  of  interpolation  errors  at  every  marching  step,  leading  to  an  apparently  divergent 
(unstable)  numerical  solution.  Increasing  the  order  of  interpolation  from  linear  to  cubic  did  not 
eliminate  the  instability. 

Looking  for  a  scheme  that  would  replicate  a  planar  CRW  accurately,  we  tried  the  modified 
marching  idea  as  presented  by  Zucrow  and  Hoffman  for  1-D  time-dependent  flows  (Sections  19-6(a) 
and  19-6(j)  of  [4] ).  In  this  scheme  new  grid  points  are  determined  by  forwardly  extending  a  "primary" 
family  of  continuous  characteristic  lines  from  old  grid  points.  The  primary  family  in  a  CRW  is  the 
characteristics  fanning  out  from  the  corner  (we  assume  it  is  the  C  family).  By  choosing  this 
modified  scheme,  the  interpolation  for  trace  points  obtained  from  reversely  extended  C  lines  was 
eliminated.  However,  the  corresponding  interpolation  for  the  transverse  C  characteristics 
remained,  and  with  it  the  aforementioned  instability. 

In  order  to  replicate  a  planar  CRW,  we  had  to  replace  the  flow  variables  by  the  Riemann 
invariants  (v±0).  In  a  C  planar  CRW,  the  Riemann  invariant  (v  +  0)  is  uniformly  constant,  so 
that  the  interpolation  in  (v  +  0)  due  to  reversely  extending  C  '  characteristics  introduces  no  error  at 
all.  This  scheme,  which  we  named  SI. VI A  (Semi  Inverse  Marching  Algorithm),  was  indeed  verified  to 
replicate  a  planar  CRW  exactly,  when  implemented  in  the  code  JET. 

The  plan  of  this  report  is  the  following.  In  Ch.  2  we  supplement  the  description  of  the  numerical 
scheme  given  in  Ch.  2  of  [1]  ,  by  adding  more  details  on  the  computational  procedure.  A  description 
of  the  code  JET  is  given  in  Ch.  3,  and  the  code  listing  is  reproduced  in  Ch.  4. 

Note  on  symmetry  : 

The  code  JET  has  two  symmetry  options.  When  DELTA  =1  a  ring-symmetry  is  in  effect;  when 
DELTA  =  0,  a  planar  symmetry  is  in  effect.  An  axisymmetric  jet  exiting  in  the  y  direction  from  the 
same  nozzle  aperture  along  the  x  axis  can  readily  be  computed  by  replacing  all  terms  in  the  code  that 
correspond  to  sin(0)/y  in  the  compatibility  equations  (2.1-1),  by  cos(0)/x.  In  that  case  the  coding  is 
virtually  unchanged,  and  the  only  care  that  should  be  exercised  is  for  the  difference  equations  for  new 
grid  points  on  or  near  the  y  axis.  Also,  all  reference  to  the  analytic  approximation  of  the  ring- 
symmetric  CRW  [1]   should  be  deleted  in  this  case,  as  it  is  designed  specifically  for  ring-symmetry. 


2.       THE  COMPUTATIONAL  SCHEME 

A  basic  description  of  the  semi  inverse  (SIMA)  and  inverse  marching  schemes  was  given  in  Ch.  2 
of  [1]  .  We  supplement  this  description  by  specifying  the  slightly  modified  definition  of  Riemann 
invariants  in  the  code,  and  by  giving  information  about  some  ancillary  computations. 


2.1      Riemann  Invariants 

The  compatibility  equations  whose  integration  constitutes  the  numerical  solution  to  the  governing 
equations  [1]   are  expressed  in  terms  of  the  Riemann  invariants  as  follows  : 

Along  C  +     (v  -  9)4  =  (v  -  6)2  +  co  sinp24  sin824  An,  /  y24 

(2.1-1) 
Along  C~     (v  +  9)4  =  (v  +  6),  +  co  sinn14  sin914  AE,  /  y,4 

The  Riemann  invariants  (v  ±  9)  are  modified  for  convenience,  by  adding  a  constant  to  both  v  and 
0.   The  new  definitions  of  v(M)  and  9  are  : 

v(M)  =    -  T  arctan(Tq)+  arctan(q) 

q  =  (M   -l)  (2.1-2) 

9  -♦  9  -  9L 

Thus,  in  a  Prandtl- Meyer  flow  with  entry  Mach  number  of  M,,  the  modified  values  of  both  v(M) 
and  9  vanish  as  M-»so.  As  a  consequence,  in  a  C  Prandtl-Meyer  flow  the  modified  invariant 
(v  +  9)  vanishes  uniformly.  In  this  modified  form,  the  computation  of  M  from  v(M)  is  readily  done 
by  performing  standard  Newton- Raphson  iterations  (in  RFUNC),  using  the  derivative  : 

v'(q)  =  -(r-oIu  +  rVXi  +  q2)!"1  (2.1-3) 


2.2     The   Integration  Scheme  for  a  New  Grid   Point 

The  integration  scheme  has  been  sketched  in  Ch.  2  of  [1]  .  It  is  performed  in  INVMAR  for  inverse 
marching  points  or  in  SEMINV  for  semi-inverse  marching  (SIMA)  points.  The  computational 
scheme  is  specified  via  the  following  seven-step  procedure  : 

INVMAR  (Inverse  Marching) 

(a)  Grid  :  At  this  stage  the  new  grid  point  has  already  been  defined. 

(b)  Predictor  :    Flow  variables  are  the  interpolated  (linear  nearest-neighbor)  value  on  the  old  row 
for  a  point  having  the  new  grid  x  coordinate  (x4). 

(c)  Centered  variables  :   Denote  the  Riemann  invariants  by 


RM  =  (v  +  6) 
RP  =  (v-6) 

then  centered  values  for  segments  (1,4)  and  (2,4)  (using  code  notation)  are  : 
RM14  =  (RMl  +  RM4)/2  RPI4  =  (RPl  +  RP4)/2 

RM24  =  (RM2  +  RM4)/2  RP24  =  (RP2  +  RP4)/2 


(2.2-1) 


(2.2-2) 


All  other  centered  flow  variables  are  computed  from  the  centered  Riemann  invariants  by  calling 
RFUNC. 

(d)  Inverse  Extension  :   old  trace  points  Xj,  X-,  are  evaluated  from  the  geometrical  relations 

Along  C"     ....     ynew  -  yoId  =  (X4-X,)  tan(914-|i14) 

(2.2-3) 
Along  C  +     ....      ynew  -  yoId  =  (x4  -  x2)  tan(624  +  n24) 

(e)  Interpolation  :   find  Riemann  invariants  RM,  RP  at  old  trace  points  Xj  and  x2  through  nearest- 
neighbor  linear  interpolation  by  calling  INTERP. 

(f)  Integration  :    Using  the  compatibility  relations  in  finite-difference  form  (2.1-1)  with  segment- 
centered  coefficients,  compute  iteration-updated  values  of  Riemann  invariants  at  new  grid  point. 


(g)      Corrector  :    if  values  of  Riemann  invariants  and  old  trace  points  X.,  ^  are  not  sufficiently 
convergent,  resume  the  procedure  at  step  (c)  above. 

SEMINV  (Semi  Inverse  Marching  -  SIMA) 

(a)  Grid  :   New  grid  point  (x4)  is  determined  as  part  of  the  SIMA  scheme  at  step  (d)  below. 

(b)  Predictor  :   Flow  variables  are  those  of  point  (x2,yo|d). 

(c)  Centered  variables  :   Identical  to  step  (c)  above. 

(d)  Semi-Inverse  Extension  :    new  grid  point  x4  and  old  trace  point  x,  are  evaluated  from  the 
geometrical  relations  in  Eq.  (2.2-3)  above. 

(e)  Interpolation  :    find  Riemann  invariants   RM,  RP  at   old  trace  point  Xj   through  nearest- 
neighbor  linear  interpolation  by  calling  INTERP. 

(f)  Integration  :   Identical  to  step  (f)  above. 

(g)  Corrector  :   Identical  to  step  (g)  above,  except  for  replacing  x2  in  the  convergence  test  by  x4. 


2.3     Boundary  Conditions 

On  the  vacuum  side  the  boundary  conditions  (p  =  0)  can  only  be  approximately  implemented  in  a 
method  of  characteristics  scheme.  We  do  so  by  ending  the  computation  on  a  certain  "final"  C  fan 
characteristic  line  that  starts  out  with  a  sufficiently  high  Mach  number  Mf  at  the  corner  (typically 
Mf=  34).  The  marching  computation  of  new  grid  points  on  the  boundary  C  characteristic  via  the 
SIMA  scheme  is  identical  to  that  of  C  '  characteristics  within  the  ring- symmetric  CRW.  It  is  noted 
that  under  this  boundary  scheme  some  outflow  takes  place  through  the  boundary  characteristic  line. 
so  that  the  total  mass  flow  throuah  a  row  y  =  y       decreases  sliahtlv  as  v       increases. 

B  j       j  new  °      •>         '  new 

At  the  nozzle  exit  the  boundarv  conditions  are  assumed  to  be  uniform  outflow  in  the  radial  (v) 
direction  with  Mach  number  M..  At  the  nozzle  lip,  the  SIMA  integration  starts  out  from  a  presumed 
planar  CRW  (Prandtl- Meyer  flow)  at  the  corner  (i.e.,  the  associate  CRW  in  the  terminology  of  Ch.  3 
in  [1]  ). 

At  the  plane  of  symmetry  (x=0)  the  boundary  condition  is  simply  9  =  7t/2.  However,  this 
condition  is  implemented  indirectly,  by  assuming  that  the  flow  at  virtual  grid  points  with  x  <  0  is  a 
mirror-image  of  the  flow  at  the  corresponding  x  >  0  points.  The  reason  is  that  when  a  new  grid 
point  of  x4=0  or  of  x4  sufficiently  close  to  zero  is  considered  for  inverse-marching  integration,  the 
inversely  extended  trace  point  (Xj,yold)  can  be  at  X  <  0.  Considering  the  subtraction  of  6L  from  0  as 
in  Eq.(2.1-2),  the  reflection  rules  are  : 


RM    -4    RP      +    (7t-2eL) 

(2.3-1) 
RP     -*    RM     -    (7t-29L) 

where  values  on  the  left  and  right  of  the  -+  symbol  correspond  to  values  left  and  right  of  x=0.   This 
boundary  condition  is  implemented  in  INTERP. 


2.4      Continuum  Breakdown  Surface 

As  an  informative  option,  the  code  JET  can  compute  (in  PLUMES)  points  on  a  surface  of 
continuum  breakdown  [5,6,7]  ,  which  is  defined  as  a  line  of  constant  B,  where  B  is  given  by  : 

B=   -(u/<p)p"l(dp/dS) 

(2.4-1) 
-1/2 

<p  =  4(7ty)       <r  n  a 

When  the  standard  isentropic  relations  for  p  and  n  in  terms  of  M  are  substituted  in  (2.4-1),  the 
flow  speed  is  expressed  as  u  =  Ma  and  the  streamwise  gradient  of  M  is  expressed  in  cartesian 
coordinates,  we  get  : 

,1/2  2    r  2,1/(7-1)-  1    r  ! 

B  =  X0(7T//8)      M    ll  +  ((Y-l)/2)M  J  IMxcosG  +Mysine] 

(2.4-2) 
^0  -  (21/2  *  njl 

Note  that  the  sign  of  B  has  been  chosen  as  positive  for  expansion  flows.  This  definition  is 
preferred  to  taking  an  absolute  value  of  the  flow  gradient,  since  it  assures  proper  interpolation  of  B 
even  if  its  spatial  distribution  goes  through  B  =  0. 

Due  to  the  dependence  of  B  on  a  spatial  gradient,  its  numerical  evaluation  (see  BREAK)  is 
attributed  to  mid-grid  points  both  in  X  and  in  y. 


3.       THE  JET  CODE 

In  this  chapter  we  provide  a  concise  description  of  the  JET  code  according  to  its  version  at  the 
time  of  the  JET018  run.  This  description  is  intended  as  an  aid  in  reading  the  code  listing  which  is 
given  in  Ch.  4. 

The  plan  of  this  chapter  is  as  follows.  Array  variables  that  constitute  the  mainstay  of  the 
computational  scheme  are  described  in  Section  3.1.  Auxiliary  array  variables  that  are  used  primarily 
for  processing  the  information  generated  by  the  numerical  scheme,  are  described  in  Section  3.2, 
followed  in  Section  3.3  by  a  list  of  major  parameters  that  control  the  computation  (some  of  them  also 
serve  as  run  data).   Finally,  all  subroutines  are  listed  and  described  in  Section  3.4. 


3.1      Main  Variables 

The  array  variables  used  for  the  computational  scheme  are  organized  in  two  labeled  COMMON 
groups.  The  first  group  /VECS/  is  designed  to  hold  two  grid  rows  -  the  old  row  designated  by  suffix 
F  and  the  new  row  designated  by  suffix  N.  The  second  group  /CHARAC/  are  characteristic-indexed 
arrays  that  hold  information  about  continuous  characteristic  lines.  This  characteristic  information  is 
used  in  two  ways  :  it  is  incorporated  in  the  SIMA  computational  scheme  for  the  CRW  region,  and  it 
is  used  to  store  data  for  optional  plotting  of  characteristic  lines  (see  PLUMES  and  PRINT). 

The  basic  organization  is  that  the  new  arrays  (suffix  N)  are  those  in  which  values  are  stored  during 
the  course  of  the  marching  computational  procedure.  At  the  end  of  each  marching  step,  values  are 
transferred  from  new  arrays  to  old  arrays  (suffix  F);  this  is  done  in  MOVE.  In  the  array  listing 
below,  we  indicate  in  parenthesis  the  subroutine  (or  subroutines)  in  which  that  new  array  is  defined. 


/VECS/ 
XN(I) 
RMN(I) 

RPN(I) 

MN(I) 

MUN(I) 

TETAN(I) 


x  coordinate  of  grid  point  I.    (GRIDN) 

modified   Riemann  invariant  (v  +  0)   at  grid  point   I.      (BEGIN,   INVMAR. 

LOADC). 

modified   Riemann  invariant  (v-6)   at  grid  point   I.      (BEGIN,    INVMAR, 

LOADC). 

Mach  number  at  grid  point  I     (BEGIN,  INVMAR,  LOADC). 

Mach  angle  \i  at  grid  point  I.     (BEGIN,  INVMAR,  LOADC). 

true  (unmodified)  flow  angle  6  at  grid  point  I.    (BEGIN,  INVMAR,  LOADC). 


BN(I)  value  of  breakdown  parameter  B  at  point  1-1/2  (and  at  half  a  marching  step 

back  in  y  as  well).    (BREAK). 
XTEMP(I)  used  for  auxiliary  computation  of  1-1/2  grid  points  in  PLUMES. 


/CHARAC/ 

XCHARN(KC) 

YCHARN(KC) 

RMCARN(KC) 

RPCARN(KC) 

TCHARN(KC) 

MUCARN(KC) 

CSIGNN(KC) 

MCHARN(KC) 
MCHARI(KC) 


X    coordinate    of  point    on    characteristic    line    number    KC.     (BEGIN, 

SEMINV,  PLUMES). 

y  coordinate  of  point  on  characteristic  line  number  KC.  (BEGIN,  SEMINV, 

PLUMES). 

modified  Riemann  invariant  (v  +  6)  of  point  on  characteristic  line  number 

KC  (BEGIN,  SEMINV). 

modified  Riemann  invariant  (v  —  0)  of  point  on  characteristic  line  number 

KC.  (BEGIN,  SEMINV). 

true  (unmodified)  flow  angle  0  at  point  on  characteristic  line  number  KC. 

(BEGIN,  SEMINV). 

Mach   angle   \i   at   point   on   characteristic   line   number   KC.      (BEGIN, 

SEMINV). 

sign  of  characteristic  line  number  KC.    It  has  value  1  for  C      and  value  —  1 

for  C  ~  .    Note  that  upon  reflection  of  a  C      line  from  the  symmetry  plane 

(X=0),  the  sign  value  is  changed  from  1  to  -  1.    (BEGIN,  SEMINV). 

Mach   number   at    point    on    characteristic    line    number    KC.      (BEGIN, 

SEMINV). 

Mach  number  at   Prandtl-Meyer's   fan  characteristic   number   KC   at   the 

corner.    It  is  defined  initially  and  is  not  changed  during  the  run.    (BEGIN). 


3.2      Auxiliary  Variables 

In  addition  to  the  major  arrays  mentioned  above,  there  are  several  groups  of  auxiliary  arrays  that 
do  not  affect  the  computational  scheme,  but  are  intended  for  informative  processing  of  the  results. 
These  groups  are  /PLUME/,  /IPLUME/,  /THICKY/,  /THICKX/.  /GRP/.  /PLUME/  is  used  to 
preserve  points  on  special  lines  for  later  plotting  (in  a  separate  code).  /THICKY/  and  /THICKX/ 
are  for  storing  values  of  radial  (y)  and  lateral  (x)  molecular  opacities.  The  group  /GRP/  is  used  in 
conjunction  with  comparative  computation  of  the  ring-symmetric  CRW  flow  according  to  the  analytic 
approximation  [1]  . 


/PLUME/      (PLUMES,  PRINT) 

XPL(J,IPL)  x  coordinate  at  marching  step  J  of  special  line  number  IPL. 

YPL(J,IPL)  y  coordinate  at  marching  step  J  of  special  line  number  IPL. 

/IPLUME/      (PLUMES,  PRINT) 

KPL  number  of  special  lines  computed  in  PLUMES. 

ITYPL(IPL)         index  indicating  the  type  of  special  line  number  IPL. 


/THICKY/      (OPACY,  PRINT) 


XTH(J) 
TH(J) 


x  coordinate  on  boundary  characteristic  line  at  marching  step  J,  from  which 
radial  opacity  is  integrated. 

radial  opacity  computed  by  y-integration  from  the  boundary  point  defined  by 
XTH(J)  (up  to  current  YN). 


/THICKX/      (OPACX,  PLUMES,  PRINT) 


YXI(JXI) 
XI(I,JXI) 

XIPM(I,JXI) 

XIGRP(I,JXI) 
XIAPP(UXI) 
XIF(I,JXI) 


y  coordinate  of  printed  row  number  JXI  (the  index  JXI  counts  just  rows  that 

have  been  printed).   The  row  to  be  printed  next  upon  calling  PRINT  is  the  row 

having  YF  near  YXI(JXI). 

lateral  (x)  molecular  opacity  [1]    at  point  XF(I),  for  printed  row  JXI.    It  is 

obtained    by    numerically    integrating    the    solution    obtained    from   the   JET 

computation  (see  OPACX). 

same  as  XI(I,JXI)  except  that  the  Prandtl- Meyer  solution  is  used  to  estimate 

the  flow  at  grid  points  XF(I). 

same  as  above,  except  that  the  analytic  approximation  to  a  ring-symmetric 

CRW  [1]   is  used  to  estimate  the  flow  at  grid  points  XF(I). 

same  as  XIGRP(I,JXI)  except  that  the  numerical  integration  is  replaced  by  an 

approximate  closed-form  expression  [1]  . 

stores  grid  points  XF(I)  of  printed  row  JXI. 


/GRP/      (PRINT,  HMSET,  MFUNC,  HINTER,  MATCH) 

DMINV  increment  of  inverse  Mach  number  for  array  MHINV(I). 

MHINV(I)  inverse   Mach  number  array  (from  0  to    1/MEXIT),  from  which  the  H(M) 

function  can  be  evaluated  (HMSET). 
HMV(I)  values  of  the  H(M)  function  evaluated  by  numerical  integration.    It  is  used  to 

compute  this  function  by  interpolation.     (HMSET,  HINTER). 


3.3      Major  Parameters 

Parameters  that  define  and  control  a  particular  run  (such  as  the  maximum  y  for  the  marching 
computation,  the  number  of  grid  points  on  a  row  and  many  more)  are  defined  in  INIDAT.  (The 
code  JET  has  no  input  file  and  no  READ  statements).  The  major  control  parameters  are  grouped  in 
/PAR/  (floating  point)  and  in  /IPAR/  (integers);  thermodynamic  data  are  grouped  in  /STAG/. 

We  indicate  in  the  listing  the  subroutines  in  which  the  labeled  COMMON  group  or  a  particular 
parameter  is  defined  (or  sometimes  referred  to). 

/PAR/      (INIDAT) 

MEXIT  nozzle  exit  Mach  number  (Mj). 

MFIN  Mach  number  of  the  final  (boundary)  CRW  characteristic  at  the  corner  (Mf). 

YiMAX  maximum  value  of  y  for  the  marching  scheme.    When  YF.GE.YMAX  the  run  is 

terminated. 

DYO  initial  marching  step. 

DY  current  marching  step. 

DYNEXT       next  marching  step  (YSTEP). 

STAB  stability  coefficient  for  marching  step  (STAB.LE.l).  (See  YSTEP). 

DELTA  symmetry  index.   DELTA  =  0  for  plane  symmetry;  DELTA  =  1  for  ring- symmetry. 

PSI1  angle  of  Prandtl- Meyer  fan  characteristic  at  exit  conditions  (measured  from  X 

axis). 

PSIF  angle  of  final  (boundary)  Prandtl-Meyer  fan  characteristic. 

SIGMA  collision  cross-section  (<r). 

FRACG  the  number  of  intervals  initially  allocated  to  the  CRW  fan  is  a  FRACG  fraction  of 

the  total  number  of  intervals  (KF0-1).    (see  BEGIN). 

EPSIL  convergence  parameter  (small  number).    (INVMAR,  SEMINV). 

TETLIM  flow  angle  (from  X  axis)  of  the  limiting  (p=  0)  velocity  vector  of  the  flow  at  the  lip- 

centered  Prandtl-Meyer  fan. 

TETSYM         PAI-2*TETLIM   for  reflection  transformation    (see  INTERP). 

/IPAR/      (INIDAT) 

JMAX  maximum  number  of  marching  steps.    If  J.GE.JMAX  run  is  terminated. 

KFO  initial  (and  maximum)  number  of  grid  points  in  a  row. 

KF  current  number  of  grid  points  in  the  old  row. 

KN  current  number  of  grid  points  in  the  new  row. 
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ITERO  maximum  number  of  iterations  for  the  integration  of  the  compatibility  relations 

(see  INVMAR  and  SEMINV;  also  used  in  RFUNC,  PLUMES). 
IM,  IP  search  indices  for  interpolation  subroutine  INTERP.    (see  INVMAR,  SEMINV). 

J  current  row  index  (also  index  of  a  marching  step). 

KF2  defined  as  2*KF;  not  used  in  present  version. 

IDEL,  JDEL  increments  for  printing  grid  point  I  and  row  J  (see  PRINT). 


JYXI 
JXI 
I  LEAD 


ILEADF 
KCLEAD 


number  of  rows  to  be  printed  in  a  run. 

index  of  printing  row,  to  be  printed  next  (see  PRINT). 

index  I  at  the  first  grid  point  on  current  new  row,  where  the  SI  MA  integration 

commences.    Initially  this  point  corresponds  to  the  leading  characteristic  of  the 

CRW.    (see  GRIDN,  BEGIN). 

value  of  I  LEAD  for  current  old  row. 

index  in  the  characteristic  array  for  the  characteristic  line  that  corresponds  to  the 

new  grid  point   I  =  ILEAD    (see  GRIDN).   Initially  KCLEAD  =  1. 


/STAG/      (INIDAT) 

RHOO,  NO       stagnation  density  and  number  density. 

PO,  TO,  AO       stagnation  pressure,  temperature  and  sound  speed. 

MDOT1  mass  flow  rate  from  ring-nozzle  (only  from  the  X  >  0  half).    (See  PRINT). 


+ 


/ICHARA/      (BEGIN) 

KCHARP        number  of  CT    characteristic  lines  for  which  data  is  stored  (either  for  SIMA 

computation  or  for  subsequent  plotting). 
KCHARM       number  of  C  ~    characteristic  lines  for  which  data  is  stored  (only  for  subsequent 

plotting). 
KCHARO        total      number      of      characteristics      for      which      data      is      stored,      i.e., 

KCHAR0=  KCHARP+  KCHARM. 
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3.4      Description  of  JET  subroutines 


MAIN  PROGRAM 

The  main  program  performs  two  functions.  The  first  section  (up  to  statement  1)  is  the  initial  set 
up;  it  is  performed  just  once.  The  second  section  is  the  marching  loop  with  the  step  index  J.  This 
program  can  be  read  as  a  flow  chart  of  the  overall  computational  procedure. 

INIDAT  is  for  setting  up  run  data.  In  BEGIN  the  initial  conditions  for  the  marching 
computation  are  set  up.  A  single  marching  step  is  performed  by  calling  MARCH,  and  the  loading  of 
new  row  vectors  into  old  row  vectors  is  done  by  calling  MOVE.  The  call  to  YSTEP  is  for  the  first 
computed  marching  step.  All  remaining  calls  are  for  informative  tasks  (see  HMSET,  BREAK, 
OPACY,  PLUMES,  PRINT).   Run  is  terminated  when  either  YF.GE.YMAX  or  when  J.GEJMAX. 

NOTE    ON    EXEC  :    The  only  special  feature  in  the  EXEC  is  retaining  the  output  unit  7  file  for 
optional  post-plotting.   The  printed  output  (unit  6)  is  the  system's  standard  (default). 


INIDAT 

Initial  data  definition  and  preliminary  data  computations.  The  data  is  defined  by  statements  rather 
than  by  reading  an  input  file.  The  meaning  of  major  parameters  was  described  in  Section  3.3  above. 
User  is  invited  to  modify  the  data  definitions,  particularly  of  run-control  parameters  such  as  YMAX, 
JYXI  and  YXI(JXI)  (for  printing  JYXI  selected  rows). 


BEGIN 


Here  all  initial  values  (prior  to  beginning  of  marching  schemes  integration)  are  loaded  into  all 
major  computational  arrays  (Section  3.1).  Also,  values  of  the  key  integer  parameters  KCHARP, 
KCHARM,  KCHARO,  ILEAD,  KCLEAD  and  KF  are  defined. 

In  the  first  loop  (loop  1)  we  define  an  initial  family  of  C  characteristic  lines  for  the  lip-centered 
CRW,    by    storing    the    Mach    number    of  the    Prandtl-Meyer    fan    characteristics    in    the    array 
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MCHARI(KC).   Note  that  the  fan  characteristics  are  generated  at  equal  RP  intervals,  since  the  flow 
variables  are  RM  and  RP.   However,  a  different  division  might  also  be  acceptable. 

The  next  step  is  the  definition  of  initial  values  for  all  characteristic  arrays,  first  the  C  arrays 
(loop  2),  then  the  C~  arrays  (loop  21).  The  C"  characteristic  lines  are  needed  just  for  informative 
output  (post-plotting),  so  the  present  version  contains  just  one  C~  line.   The  user  may  modify  that. 

The  remaining  grid  points  (altogether  KFO  grid  points  are  initially  available)  are  uniformly 
distributed  across  the  nozzle  opening,  and  the  row  arrays  are  loaded  with  the  corresponding  nozzle- 
exit  flow  variables  (loop  3). 


PRINT 

The  main  task,  of  this  subroutine  is  the  printing  of  flow  variables  at  grid  points  of  selected  rows. 
The  printing  of  a  row  is  selected  when  YF  is  close  to  a  predefined  array  YXI(JXI).  Following  the 
printing,  JXI  is  updated  by  adding  1. 

For  comparison,  additional  flow  variables  are  printed  for  each  row.  These  are  computed  from  the 
analytic  approximation  to  a  ring- symmetric  CRW  [1]  ,  by  calling  MATCH.  Also,  lateral  molecular 
opacities  of  various  kinds  of  approximation  are  computed  by  calling  OPACX,  and  are  printed  for 
each  grid  point  within  the  CRW. 

Following  the  row  printing  (statement  120),  arrays  intended  for  post-processing  (plotting  of  special 
lines)  are  printed  and  subsequently  written  on  output  unit  7.  This  is  done  once  per  run,  just  before 
run  termination. 


FIN 


This  subroutine  is  called  when  an  error  is  encountered,  in  order  to  terminate  the  run.  Note  that 
the  run  is  terminated  by  deliberately  introducing  an  error  of  computing  SQRT(-l),  which  is  done  in 
order  to  trigger  the  printing  of  calling  sequences  by  the  operating  system. 
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MARCH 

This  subroutine  performs  a  single  marching  step  by  calling  the  proper  computational  subroutines 
at  an  appropriate  sequence.  It  can  be  read  as  a  flow  chart  of  the  entire  computational  scheme.  First 
the  segment  of  the  new  row  suitable  for  SI  MA  computation  is  calculated  by  calling  SEMINV.  Then 
new  grid  points  for  that  part  of  the  new  row  for  which  inverse  marching  integration  is  to  be 
performed,  are  generated  by  calling  GRIDN.  The  results  of  the  SEMINV  computation,  which  were 
stored  in  characteristic  arrays,  are  now  loaded  into  row  arrays  by  calling  LOADC.  Finally,  the 
computation  of  the  new  row  is  completed  by  calling  INVMAR  which  computes  the  flow  at  the 
remaining  grid  points  by  the  inverse  marching  scheme. 


INVMAR 

This  is  one  of  the  two  central  subroutines  for  computing  the  flow  at  new  grid  points  (the  other  is 
SEMINV).  Here  the  inverse  marching  scheme  is  used.  The  computational  procedure  follows  the 
seven-step  description  given  in  Section  2.2  above.  Note  that  the  initial  value  of  the  search  indices  IM 
and  IP  is  not  redefined  at  each  call  to  INTERP,  since  it  is  assumed  that  IM  and  IP  do  not  change 
much  at  consecutive  calls  to  INTERP,  so  that  search  efficiency  is  enhanced  by  not  starting  the  search 
from  an  arbitrary  point  (such  as  either  end  of  the  row). 


SEMINV 

This  is  the  subroutine  performing  the  SI  MA  scheme  for  computing  the  flow  at  new  grid  points 
located  along  continuous  characteristic  lines  of  the  lip-centered  CRW  (at  prescribed  y-marching 
steps).  The  essence  of  the  computational  procedure  of  this  subroutine  was  given  as  a  seven-step 
description  in  Section  2.2.  The  same  remark  about  IM  given  in  the  preceding  INVMAR  description 
applies  here  as  well. 

The  main  loop  (100)  is  over  all  characteristic  lines,  including  some  C  "  lines  in  addition  to  the  C 
lines.   Thus,  the  array  CSIGNF(KC)  is  used  to  get  the  appropriate  expressions  for  either  C      or  C~ 
characteristics.    It  is  noted  that  while  normally  the  characteristic  segments  through  points  1  and  2  are 
C      and  C      respectively,  this  is  reversed  when  aC"    rather  than  a  C      line  is  computed  via  the 
SIMA  scheme.    In  this  case,  which  is  characterized  by  having  CSIGNF(KC).LT.O,  the  Riemann 
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invariants  integrated  along  segments  (1,4)  and  (2,4)  are  interchanged.    This  is  done  in  the  few 
statements  just  preceding  and  following  statement  21. 

An  additional  capability  of  this  subroutine  is  to  treat  a  change  of  a  C  characteristic  line  into  a 
C  ~~  line  upon  reflection  from  the  symmetry  plane  (x  =  0).  This  is  done  by  first  computing  a  new  grid 
point  having  X4.LT.0,  and  then  changing  its  sign  after  setting  CSIGNN(KC)=  -  1  (statements  just 
preceding  statement  30).  It  is  also  possible  to  skip  the  computation  of  a  particular  characteristic  by 
setting  CSIGNN(KC)  =  0.   This  feature  is  not  exploited  in  the  present  version. 

Finally,  we  note  that  not  all  characteristic  lines  computed  here  are  part  of  the  marching  flow 
computation.  Only  those  with  indices  KC  between  KCLEAD  and  KCHARP  are.  All  other 
characteristic  lines  are  computed  just  for  informative  purposes  (post-plotting). 


RFUNC 

Here  M,  MU,  TETA  are  computed  from  the  two  Riemann  invariants  RM,  RP.  The  computation 
of  M  is  performed  by  a  Newton-Raphson  iteration  using  Equations  (2.1-2)  and  (2.1-3)  given  in 
Section  2.1  above. 


INTERP 

This  subroutine  starts  by  finding  through  a  search  procedure  the  grid  interval  (I,  1  +  1)  that 
contains  a  given  point  X.  Then  the  Riemann  invariants  are  computed  for  this  point  by  linear 
interpolation,  and  returned  in  RM,  RP.  Note  that  X  may  be  negative,  which  accounts  for  the 
relatively  elaborate  search  logic  in  the  determination  of  I,  and  for  the  reflection  transformation  (as  in 
Eq.  (2.3-1)  above)  preceding  the  last  two  statements  of  the  subroutine. 


INTERX 

This  interpolation  routine  performs  an  inverse  task  to  that  of  INTERP,  in  that  it  finds  the  point 
X0  that  corresponds  to  a  given  linearly  interpolated  value  of  the  flow  variable  VARO.  It  is  used  in 
PLUMES  to  compute  the  location  of  a  breakdown  surface  point  on  a  new  row  of  x-centered  and  y- 
centered  grid  points 
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BREAK 

This  subroutine  computes  the  new  breakdown  parameter  array  BN(I).   The  computation  is  based 
on  the  description  given  in  Section  2.4  above. 


OPACY 

Here  the  radial  (Y)  molecular  opacity  array  TH(J)  is  computed.  At  each  marching  step  J,  a  new 
boundary  grid  point  XTH(J)  is  added,  then  the  radial  opacities  at  all  preceding  boundary  points  are 
updated  by  adding  the  contribution  of  the  gas  layer  between  the  current  old  and  new  rows.  Note  that 
since  grid  points  on  adjacent  rows  are  not  located  on  equal-x  columns,  this  procedure  requires  x- 
interpolation  by  calling  INTERP. 


PLUMES 

This  is  a  user-defined  subroutine,  where  up  to  10  special  lines  can  be  computed  and  subsequently 
retained  on  output  unit  7  for  post-processing  (plotting).  The  type  of  the  line  ITYPL(IPL)  and  a 
parameter  VPL(IPL)  that  defines  the  line,  are  computed  through  user-inserted  statements  in  the 
section  preceding  statement  2000.  Then  an  additional  point  on  the  current  new  row  is  computed  for 
each  line  type.  The  available  types  are  clearly  stated  in  comments.  Note  that  characteristic  lines 
have  already  been  computed  in  SEMINV  using  the  SI  MA  scheme,  regardless  of  whether  they  are  part 
of  the  solution  grid  to  the  flow  field,  or  are  just  computed  for  informative  purpose.  It  is  the  user's 
choice  which  of  these  lines  (if  any)  are  to  be  saved  in  the  /PLUME/  arrays  for  subsequent  post- 
processing (plotting). 


GRIDN 

This  subroutine  computes  the  grid  points  in  that  segment  of  the  new  row  for  which  the  flow  is 
computed  by  the  inverse  marching  scheme  (in  INVMAR).  Initially,  this  segment  extends  from  x=0 
to  the  new  row  grid  point  which  lies  on  the  leading  characteristic  of  the  lip-centered  CRW.  However, 
since  the  leading  characteristic  is  reflected  from  the  symmetry  plane  (x  =  0)  at  some  point,  this 
segment  steadily  shrinks  in  size  as  the  marching  proceeds.    The  remedy  is  to  declare  the  next-to-the- 
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leading  characteristic  line  (KC  =  2)  as  the  beginning  of  the  segment  for  SI  MA  integration,  by  setting 
KCLEAD  =  2.  This  process  of  increasing  KCLEAD  is  repeated  whenever  it  is  deemed  necessary. 
The  criterion  in  the  present  version  for  the  minimal  KCLEAD  is  that  the  inverse-marching  segment 
should  be  at  least  twice  DX1  -  the  average  CRW  grid  interval  (loop  1,  the  two  statements  following 
DX1  =  ...).  Also,  I  LEAD  is  redefined  for  each  row  according  to  XLEAD/DX1  +  2  in  order  to  achieve 
a  row  of  relatively  uniform  grid  intervals  throughout.  The  result  is  that  the  number  of  grid  points  in 
a  row  is  initially  KFO,  but  eventually  it  decreases  due  to  both  increase  of  KCLEAD  and  decrease  of 
I  LEAD. 


YSTEP 

In  this  subroutine  the  next  marching  step  DYNEXT  is  computed  at  the  end  of  the  current 
marching  step.    It  is  defined  as  the  smallest  step  obtained  by  forward  intersection  of  C      and  C 
characteristics  from  adjacent  grid  points.    Note  that  the  actual  value  of  DYNEXT  is  reduced  by  a 
"stability"  factor  STAB,  and  that  DY  is  also  limited  by  the  growth-rate  factor  DDY  and  by  DYMAX 
(see  MAIN  PROGRAM). 


MOVE 

Here  old  row  arrays  (loop  1)  and  old  characteristic  arrays  (loop  2)  are  loaded  with  values  of  flow 
variables  from  corresponding  new  arrays,  in  preparation  for  the  next  marching  step.  As  a  result  of 
this  organizational  feature,  informative  computations  (e.g.  BREAK,  OPACY)  that  require  both  new 
and  old  rows,  have  to  be  performed  prior  to  calling  MOVE. 


OPACX 

Here  lateral  (X)  opacities  that  correspond  to  the  number  of  expected  collisions  of  a  fast  molecule 
invading  the  CRW  in  the  —X  direction,  are  computed.  All  opacities,  except  XIAPP(I),  are  computed 
by  numerical  integration.  In  loop  1  we  compute  the  opacity  contribution  of  the  segment  lying  just 
outside  the  computational  boundary  characteristic  (MFIN),  assuming  a  Prandtl-Meyer  flow.  This 
additional  opacity  is  denoted  XIO.  If  the  flow  is  ring-symmetric,  XIO  is  recalculated  using  the  analytic 
approximation  [1]  to  estimate  the  flow  field  at  the  fringes  of  the  ring- symmetric  CRW  (see  also  the 
closed  form  expression  for  T  in  [1]  ). 
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The  computation  of  opacity  arrays  starts  after  statement  14.  First,  the  opacity  at  each  grid  point 
is  set  to  XIO.  Thus,  even  though  the  numerical  flow  computation  does  not  include  the  fluid  outside 
the  boundary  characteristic  line,  the  opacity  integration  includes  an  estimate  of  that  "missing"  part, 
i.e.,  of  XIO.  In  typical  case  computations  of  a  ring- symmetric  CRW  [1]  we  found  that  the  maximum 
value  of  XIO  was  about  0.16.,  which  indicated  that  as  far  as  interaction  with  invading  ambient 
molecules  is  concerned,  the  approximation  MFIN=  34  was  a  reasonable  substitute  for  MFIN=  oo. 

The  next  step  is  the  computation  by  numerical  integration  of  three  approximations  to  the  lateral 
opacity  :  XI(I,JXI),  XIPM(I,JXI),  XIGRP(I,JXI).  (Note  that  when  the  flow  is  ring-symmetric,  the 
approximation  XIPM(I,JXI)  obtained  by  assuming  a  Prandtl-Meyer  flow  is  usually  grossly 
exaggerated).  The  opacity  XIGRP(I,JXI)  is  based  on  the  analytic  approximation  to  a  ring- symmetric 
CRW  [1]  ,  and  is  reasonably  close  to  XI(I,JXI)  which  is  obtained  from  the  numerical  solution  to  the 
flow  field.  Finally,  a  simplified  closed-form  integration  of  lateral  opacity  [1]  is  computed  as 
XIAPP(UXI)  (loop  3).  Thus,  the  quantitative  difference  between  XI(I,JXI)  and  XIGRP(I,JXI)  is  an 
indication  to  the  degree  of  accuracy  achieved  by  the  analytic  approximation  to  a  ring- symmetric 
CRW  [1]  ,  while  the  difference  between  XIGRP(I,JXI)  and  XIAPP(I,JXI)  indicates  the  level  of  error 
introduced  by  the  closed-form  integration  of  lateral  opacity  [1]  . 


LOADC 


Here  flow  variables  of  new  grid  points  computed  via  the  SIMA  scheme  (SEMINV)  are  loaded  into 
new  row  arrays  from  corresponding  characteristic  arrays. 


NUFUNC 

This  function  computes  the  modified  v(M)  value  as  given  by  Eq.  (2.1-2).    Note  that  presently 
NU0  =  0(see  INTDAT). 


HMSET 

This  subroutine  is  called  just  once  from  the  MAIN  PROGRAM.  Its  task  is  to  set  up  the  arrays  in 
/GRP/,  so  that  the  function  H(M)  [1]  can  be  evaluated  by  interpolation  (in  HINTER).  There  is  also 
an  informative  printout  of  various  derivatives  (see  Ch.  3  of  [1]  )  generated  in  this  subroutine. 
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MFUNC 

This  subroutine  is  called  by  HMSET  in  order  to  compute  functions  of  Mach  number  that  serve  in 
the  computation  of  H(M).  The  output  variable  F  is  the  integrand  for  the  integration  leading  to 
H(M). 


HINTER 

This  subroutine  computes  H(M)  by  linear  interpolation  in  inverse  Mach  number,  using  the  /GRP/ 
arrays  computed  in  HMSET. 


MATCH 

This  subroutine  is  called  from  PRINT  to  compute  the  Mach  number  according  to  the  analytic 
approximation  of  a  ring- symmetric  CRW  [1J  ,  for  point  (YF,XF(I)).  MOB  is  the  associate  Mach 
number  M(0,p),  which  is  preserved  in  the  array  MCHARI(KC)  for  all  CRW  characteristics  that  are 
used  in  the  SIMA  computation.  Hence  the  Mach  number  M(a,p),  denoted  by  MAB  can  be 
computed  directly  from  the  analytic  approximation  [1]  to  the  area  function  at  (YF,XF(I))  by  calling 
AREAF.  Since  typically  M(O.P)  is  not  known,  we  also  compute  the  Mach  number  via  the  inverse- 
problem  procedure  [1]  ,  denoting  the  resulting  Mach  numbers  by  suffix  I  :  MOBI  for  M(O.P)  and 
MAB  I  for  M(a,P).  The  inverse-problem  iterative  procedure  [1]  is  performed  in  loop  1.  resulting  in 
MOBI.  From  MOBI  the  value  of  MABI  is  computed  through  the  area  function  approximation  as  for 
MAB  above. 


AREAF 

This  subroutine  computes  the  Mach  number  M  that  corresponds  to  the  area  function  F  (Eq. 
(3.2-1)  of  [1]  ).  The  computation  is  done  by  Newton- Raphson  iterations,  and  it  has  been  found  to 
converge  when  M.GT.l  (and  when  M  —  1  is  not  much  smaller  than  1). 
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4.       THE  JET  CODE  LISTING 


CSOPTIONS  LIST 

C  JET018 

C   "JET"   A  SEMI-INVERSE  MARCHING  CHARACTERISTICS  METHOD  FOR  RING  JETS. 

C   USING  RIEMANN  INVARIANTS  RM=(NU+TETA) ,  RP=(NU-TETA)  AS  FIELD 

C   VARIABLES. 

IMPLICIT  REAL*8(A-H,L-Z,$) 

REAL**  XPL,YPL 

COMMON  /PLUME/XPL ( 1002, 10 ),YPL( 1002) 

COMMON  /IPLUME/KPL,ITYPL(10) 

COMMON  /VECS/XF(101),RMF(101),RPF(101),MF(101),MUF(101), 

1  TETAF(101),BF(101), 

2  XN(101),RMN(101),RPN(101),MNa01),MUN(101), 

3  TETAN(101),BN(101),XTEMP(101) 
COMMON/THICKY/XTHC1 002), TH( 1002) 

REALMS  YXI,XI,XIPM,XIGRP,XIAPP,XIF 

COMMON  /THICKX/YXI(20),XI(101,20),XIPM(101,20),XIGRP(101,20) 
1  ,XIAPP(101,20),XIFC101,20) 

COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15 
1  G16,G17,G18,G19,G20 

COMMON  /PAR/PAI,PAI2,DEG,XC,YC,MEXIT,MFIN,YMAX,DY0,DY,DYNEXT, 

1  STAB, DELTA, PSI1 , PSIF,ZETA1 , SIGMA, FRACG, EPSIL,NU0, 

2  TETSYM,TETLIM,DDY,DYMAX 
COMMON  /STAG/RHO0,N0,P0,T0,A0,MDOTl 
COMMON  /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP,J, 

1  KF2,IDEL,JDEL,JYXI,JXI,ILEAD,ILEADF,KCLEAD 

COMMON  /ROW/YF,YN,DXF,DXN 
COMMON  /CHARAC/XCHARF(92),YCHARF(92),XCHARN(92),YCHARN(92), 

1  RMCARFC  92) , RPCARFC  92) , RMCARNC92) , RPCARNC92)  , 

2  TCHARF(92),TCHARN(92),MUCARF(92),MUCARN(92), 

3  CSIGNN(92),CSIGNF(92),MCHARN(92),MCHARF(92), 

4  MCHARK92) 

COMMON  /ICHARA/KCHARP, KCHARM, KCHARO 
COMMON  /GRP/DMINV,MHINV(101),HMV(101) 
COMMON  /IGRP/KHM 
C 

PRINT  101 
101   FORMATCl') 
J  =  l 

IF(J.EQ.l)  STOP 
CALL  INIDAT 
PRINT  101 
CALL  HMSET 
PRINT  101 
CALL  BEGIN 

MARCH 

OPACY 

PLUMES 

PRINT 


CALL  TO  GRIDN 


CALL 

CALL 

CALL 

CALL 

J=2 

CALL  PLUMES 

CALL  MOVE 

CALL  OPACY 

CALL  PRINT 

CALL  YSTEP 

J  =  J  +  1 
DY  WAS  DETERMINED  BY  THE  PREVIOUS 

DY=DMIN1( DYNEXT, DY*DDY, DYMAX) 
INTEGRATE  BY  ONE  Y-STEP 

CALL  MARCH 
BREAKDOWN  PARAMETER  (BF(D). 

CALL  BREAK 
SPECIALLY  DESIGNATED  LINES  (FOR  PLOTTING) 

CALL  PLUMES 
STORE  NEW  LINE  (N)  IN  OLD  LINE  (F). 

CALL  MOVE 
COMPUTE  RADIAL 

CALL  OPACY 
Y-STEP  IS  VARIABLE,  SO  JMAX 

IF(YF.GE.YMAX)  JMAX=J 
PRINT  FIELD  AT  MOST  RECENT  Y 

CALL  PRINT 
NEXT  Y-STEP. 


MOLECULAR  OPACITIES 


IS  USED  AS  END-OF-RUN  CRITERION. 


JET0001 
JET0002 
JET0003 
JET0004 
JET0005 
JET0006 
JET0007 
JET0008 
JET0009 
JET0010 
JET0011 
JET0012 
JET0013 
JET0014 
JET0015 
JET0016 
JET0017 
,JET0018 
JET0019 
JET0020 
JET0021 
JET0022 
JET0023 
JET0024 
JET0025 
JET0026 
JET0027 
JET0028 
JET0029 
JET0030 
JET0031 
JET0032 
JET0033 
JET0034 
JET0035 
JET0036 
JET0037 
JET0038 
JET0039 
JET0040 
JET0041 
JET0042 
JET0043 
JET0044 
JET0045 
JET0046 
JET0047 
JET0048 
JET0049 
JET0050 
JET0051 
JET0052 
JET0053 
JET0054 
JET0055 
JET0056 
JET0057 
JET0058 
JET0059 
JET0060 
JET0061 
JET0062 
JET0063 
JET0064 
JET0065 
JET0066 
JET0067 
JET0068 
JET0069 
JET0070 
JET0071 
JET0072 
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:  JETPR 


FORTRAN   Al 


CALL  YSTEP 

IF(J.LT.JMAX)  GO  TO  1 
STOP 
FNH 


IHU>A-T^ 


ET0073 
ET0074 
ET0075 

ET0076 


SUBROUTINE  INIDAT  JET0077 

SUBROUTINE  NUMBER   1  JET0078 

IMPLICIT  REAL*8CA-H,L-Z,$)  JET0079 

COMMON  /VECS/XFC101),RMF(101),RPF(101),MF(101),MUF(101),  JET008  0 

1  TETAF(101),BF(101),  JET0081 

2  XN(101),RMN(101),RPN(101),MN(101),MUN(101),  JET0082 

3  TETAN(101),BN(101),XTEMPU01)  JET0083 
COMMON/THICKY/XTHC1 002),  THC  1002)  JET0084 
REAL**  YXI,XI,XIPM,XIGRP,XIAPP,XIF  JET0085 
COMMON  /THICKX/YXI(20),XI(101,20),XIPM(101,20),XIGRP(101,20)  JET0086 

1  ,XIAPP(101,20),XIF(101,20)  JET0087 
COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15,JET0088 

1              G16,G17,G18,G19,G20  JET0089 

COMMON  /PAR/PAI,PAI2,DEG,XC,YC,MEXIT,MFIN,YMAX,DY0,DY,DYNEXT,  JET0090 

1  STAB, DELTA, PSI1 , PSIF, ZETA1 , SIGMA, FRACG, EPSIL , NUO,  JET0091 

2  TETSYM,TETLIM,DDY,DYMAX  JET0092 
COMMON  /STAG/RHO0,N0,P0,T0,A0,MDOTl  JET0093 
COMMON  /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP,J,  JET0094 

1              KF2,IDEL,JDEL,JYXI,JXI,ILEAD,ILEADF,KCLEAD  JET0095 

COMMON  /ROW/YF,YN,DXF,DXN  JET0096 

JET0097 

PAI=4.D0*DATANC1.D0)  JET0098 

PAI2=2.D0*DATAN(1.D0)  JET0099 

DEG=180.D0/PAI  JET0100 

AR=8.3143D3  JETO101 

AV=6.022D  26  JET0102 

AW=7.27D0  JET0103 

RHO0=0.0075D0  JET0104 

T0=2300.D0  JET0105 

G=1.54D0  JET0106 

D=2.5D-10  JET0107 

MEXIT=4.D0  JET0108 

MFIN=34.D0  JET0109 

XC=0.5D0  JET0110 

YC=2.5D0  JET0111 

DELTA=0  CORRESPONDS  TO  PLANE  SYMMETRY  JET0112 

DELTA=1  CORRESPONDS  TO  CYLINDRICAL  SYMMETRY  JET0113 

DELTA=1.D0  JET0114 

FRACG=0.6D0  JET0115 

EPSIL=l.D-8  JET0116 

ITER0=20  JET0117 

KF0=101  JET0118 

JMAX=1001  JET0119 

STAB=0.50D0  JET0120 

DDY=1.05D0  JET0121 

DYMAX=0.5D0  JET0122 

YMAX=50.D0  JET0123 

DY0=YC/250.D0  JET0124 

IDEL=1  JET0125 

JDEL=1  JET0126 

POINTS  FOR  PRINTING  FLOW  FIELD  AT  YF=YXI(JXI)  JET0127 

JXI=1  JET0128 

JYXI=11  JET0129 

DYXI=5.D0  JET0130 

YXI(1)=YC+0.5D0  JET0131 

YXI(2)=YXI(1)+2.D0  JET0132 

10=2  JET0133 

DO  1  1=10, JYXI  JET0134 

YXI(I)=YXI(I0)+DYXI*DFLOAT(I-I0)  JET0135 

CONTINUE  JET0136 

IFCKF0.GT.101)  CALL  FIN(lOl)  JET0137 

IF(JMAX.GT.lOOl)  CALL  FINC102)  JET0138 

IF(FRACG.GT.1.D0  .OR.  FRACG. LT.O.)  CALL  FINC103)  JET0139 

IFCJYXI.GT.20)  CALL  FINC104)  JET0140 

IF(DELTA*(1. DO-DELTA) .NE.O.)  CALL  FIN(105)  JET0141 

N0=RHO0*AV/AW  JET0142 

AO=DSQRT(G*AR*TO/AW)  JET0143 

P0=AR*RHO0*T0/AW  JET0144 
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JETPR     FORTRAN   Al 

SIGMA=PAI*D*X2  JET0145 

LAMDA0=1.D0/(DSQRT(2.D0)XSIGMA*N0)  JET0146 

Gl=(G-l.DO)/2.DO  JET0147 

G2=(G+1.D0)/(2.D0XCG-1.D0))  JET0148 

G3=G/2.DO  JET0149 

G4=(G+1.D0)/(G-1.D0)  JET0150 

G5=DSQRT((G+1.D0)/(G-1.D0))  JET0151 

G6=1.D0/(G-1.D0)  JET0152 

G7=2.D0/(G+1.D0)  JET0153 

G8  =  (0.5D0X(G+1.D0)XX2/(G-1.D0))XX(1„D0/(G+1.D0))X  JET0154 

1    (CG+1.D0V(G-1.D0))XX((G-1.D0)/(G+1.D0))  JET0155 

G9=(G+3.D0)/C2.D0X(G-1.D0))  JET0156 

GlO=(7.DO-3„D0XG)/(2.DOX(G-l.DO))  JET0157 

G11=(2.D0/(G+1.D0))XX(1.D0/CG-1.D0))  JET0158 

G12=DSQRT((G+l.DO)/(G-l.DO))-l.DO  JET0159 

G13=(2.DO-G)/(2.DOX(G-l.DO))  JET016  0 

G14=G/(2.D0X(G-1.D0))  JET0161 

G15=(G+1.D0)/(3.D0-G)  JET0162 

G16=(G+1.D0)/4.D0  JET0163 

G20=LAMDA0XDSQRT(PAIXG/8.D0)  JET0164 

ZETA1=G5XDATAN(DSQRT(MEXITXX2-1 .DO)/G5)  JET0165 

AMUl=DARSIN(l.DO/MEXIT)  JET0166 

PSI1=PAI2+AMU1  JET0167 

ZETAF=G5XDATAN(DSQRT(MFINXX2-l.DO)/G5)  JET0168 

PSIF=PSI1+ZETA1-ZETAF  JET0169 

NUO=0.  JET0170 

TETLIM=NUFUNC(MEXIT)+PAI2-NUO  JET0171 

PSILIM=TETLIM  JET0172 

TETSYM=PAI-2.DOXTETLIM  JET0173 

GOREM=1.DO+G1XMEXITXX2  JET0174 

RH01=RHOO/GOREMXXG6  JET0175 

V1=MEXITXAO/DSQRT(GOREM)  JET0176 

P1=PO/GOREMXX(G/(G-1.DO))  JET0177 

Tl=T0/DSQRT(GOREM)  JET0178 

YYC=2.D0XPAIXYC  JET0179 

IFCDELTA.EQ.O.)  YYC=l.DO  JET0180 

MD0T1=YYCXXCXRH01XV1  JET0181 

C  JET0182 

PRINT  21,AR,AV,AW,G,RHO0,N0,P0,T0,A0,D  JET0183 

21  F0RMAT(/1X, 'THERMODYNAMIC  DATA: •/  JET0184 

1  IX, »AR,AV,AW,G=',2X,2D14.5,2F9.3/  JET0185 

2  IX, 'RHO0,N0,P0,T0,A0,D=',6D13.5)  JET0186 
PRINT  22,XC,YC,MEXIT,RH01,P1,T1,V1,MD0T1,PSI1XDEG,PSIFXDEG,  JET0187 

1          PSILIMXDEG  JET0188 

22  FORMATC/1X, 'CORNER  DATA:   XC, YC= ' , 2F9 . 2/  JET0189 

1  IX, 'EXIT  CONDITIONS: ',  JET0190 

2  2X,»MEXIT,RH01,P1,T1,V1,MD0T1  =  SF9.3,5D13.4/  JET0191 

3  1X,'CENTERED  FAN  LIMITS:',  JET0192 

4  2X, 'PSI1,PSIF,PSILIM=',3F10.3)  JET0193 
PRINT  23, DELTA, KFO , JMAX, ITERO, DYO, YMAX, STAB, DDY  JET0194 

23  FORMATC/1X, 'INTEGRATION  DATA.       SYMMETRY  INDEX:   DELTA= ' , F4 . 1/  JET0195 

1  IX, 'NUMBER  OF  POINTS  IN  X  AND  Y  DIRECTIONS:  KFO,JMAX=',   JET0196 

2  215/  JET0197 

3  IX, 'MAX.  NUM.  OF  ITERATIONS    ITERO=»,I5/  JET0198 

5  IX, 'INITIAL  Y-STEP  AND  MAXIMUM  Y:   DYO , YMAX= ' , 2D14 . 5/  JET0199 

6  IX, 'Y-STEP  STABILITY  FACTORS  STAB, DDY= • , 2F7 . 3)  JET0200 
RETURN  JET0201 
END                                                 &£GlN  JET0202 

SUbKUUIlNh  BEGIN JEIU2U3 

C   SUBROUTINE  NUMBER   2  JET0204 

IMPLICIT  REAL*8(A-H,L-Z,$)  JET0205 

COMMON  /VECS/XF(101),RMFC101),RPF(101),MF(101),MUF(101),  JET0206 

1  TETAF(101),BF(101),  JET0207 

2  XN(101),RMNU01),RPN(101),MN(101),MUN(101),  JET0208 

3  TETAN(101),BN(101),XTEMP(101)  JET0209 
COMMON/THICKY/XTH(l 002), TH( 1002)  JET0210 
COMMON  /GAMA/G,G1,G2,G3,G<+,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15,JET0211 

1              G16,G17,G18,G19,G20  JET0212 

COMMON  /PAR/PAI,PAI2,DEG,XC,YC,MEXIT,MFIN, YMAX, DYO, DY,DYNEXT,  JET0213 

1  STAB, DELTA, PSI1 , PSIF, ZETA1 , SIGMA, FRACG, EPSIL , NUO,  JET0214 

2  TETSYM,TETLIM,DDY,DYMAX  JET0215 
COMMON  /STAG/RHOO,NO,PO,TO,AO,MDOT1  JET0216 
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C 
C 
C 


21 


COMMON  /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP,J,  JET0217 

1             KF2,IDEL,JDEL,JYXI,JXI,ILEAD,ILEADF,KCLEAD  JET0218 

COMMON  /ROW/YF,YN,DXF,DXN  JET0219 

COMMON  /CHARAC/XCHARFC92),YCHARFC92),XCHARNC92),YCHARNC92),  JET0220 

1  RMCARFC92),RPCARFC92),RMCARNC92),RPCARNC92),  JET0221 

2  TCHARFC92),TCHARNC92),MUCARFC92),MUCARNC92),  JET0222 

3  CSIGNNC92),CSIGNFC92),MCHARNC92),MCHARFC92),  JET0223 

4  MCHARK92)  JET0224 
COMMON  /ICHARA/KCHARP,KCHARM,KCHARO  JET0225 

JET0226 

DEFINE  INITIAL  CHARACTERISTIC  PARAMETERS.  USE  INTERPOLATION  OF       JET0227 

RIEMANN  INVARIANT  ACROSS  THE  FAN.  JET0228 

KCHARP=IDINTCFRACG*DFL0ATCKF0-l)+l.D-6)+l  JET0229 

KCHAR0=KCHARP+1  JET0230 

KCHARM=KCHARO-KCHARP  JET0231 

IFCKCHARP.LT. 2  )  CALL  FINC200)  JET0232 

IFCKCHAR0.GT.92)  CALL  FINC210)  JET0233 

IFCKCHARM.LT.  1)  CALL  FINC205)  JET0234 

NU1=NUFUNCCMEXIT)  JET0235 

RM1=NU0  JET0236 

TET1=RM1-NU1  JET0237 

RP0=NU1-TET1  JET0238 

NUFIN=NUFUNCCMFIN)  JET0239 

RPFIN=NUFIN-CRM1-NUFIN)  JET0240 

DRP=CRPFIN-RP0)/DFLOATCKCHARP-l)  JET0241 

DO  1  KC=1,KCHARP  JET0242 

RPl=RPO+DRP*DFLOATCKC-l)  JET0243 

CALL  RFUNCCRM1,RP1,M1,MU1,TETA1)  JET0244 

MCHARICKC)=M1  JET0245 

CONTINUE  JET0246 

DATA  FOR  C+  CHARACTERISTICS.  JET0247 

THE  RIEMANN  INVARIANTS  ARE  DEFINED  IN  SUCH  A  WAY  THAT  BOTH  VANISH  AT  JET0248 

INFINITE  MACH  NUMBER.  JET0249 

RM1=NU0  JET0250 

DO  2  KC=1,KCHARP  JET0251 

CSIGNFCKC)=1.D0  JET0252 

XCHARFCKC)=XC  JET0253 

YCHARFCKC)=YC  JET0254 

IFCMCHARKKO.EQ.O.)  CALL  FINC231)  JET0255 

NU  =  NUFUNCCMCHARICKO)  JET0256 

TET=RM1-NU  JET0257 

RP1=NU-TET  JET0258 

CALL  RFUNC(RM1,RP1,M1,MU1,TETA1)  JET0259 

MCHARFCKC)=M1  JET0260 

MUCARFCKC)=MU1  JET0261 

TCHARFCKC)=TETA1  JET0262 

RMCARFCKC)=RM1  JET0263 

RPCARFCKC)=RP1  JET0264 

CONTINUE  JET0265 

DATA  FOR  C-  CHARACTERISTICS.  JET0266 

KC1=KCHARP+1  JET0267 

XCHARFCKC1)=0.8D0*XC  JET0268 

DO  21  KC=KC1,KCHAR0  JET0269 

CSIGNF(KC)=-1.D0  JET0270 

MCHARICKC)=MEXIT  JET0271 

MUCARFCKC)  =  DARSINC1.DO/MCHARICKO)  JET027  2 

TCHARFCKC)=PAI2  JET0273 

YCHARFCKC)=YC  JET0274 

MCHARFCKC)=MEXIT  JET0275 

RMCARFCKC)=RM1  JET0276 

RPCARFCKC)=NUFUNCCMEXIT)-CTCHARFCKC)-TETLIM)  JET0277 

CONTINUE  JET0278 

DEFINE  GRID  AND  INITIAL  CONDITIONS  AT  EXIT  PLANE.                    JET0279 

KFAN=KCHARP-1  JET0280 

ILEAD=KFO-KFAN  JET0281 

KCLEAD=1  JET0282 

KF=KF0  JET0283 

KF2  =  2*KF  JET0284 

YF=YC  JET0285 

DO  3  1=1, KF  JET0286 

KC=KCLEAD+I-ILEAD  JET0287 

IFCKC.GT.KCHARP)  CALL  FINC241)  JET0288 
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IFCKC.GE.l)  GO  TO  31 

XF(I)=DFLOAT(I-l)*XC/DFLOATCILEAD-l) 

MF(I)=MEXIT 

TETAF(I)=PAI2 

GO  TO  32 

31  CONTINUE 
XF( I )=XC 
MF(I)=MCHARF(KC) 
TETAF(I)=TCHARF(KC) 

32  CONTINUE 
RMF(I)=NUFUNC(MF(I))+(TETAF(I)-TETLIM) 
RPF(I)=NUFUNC(MF(I))-(TETAF(I)-TETLIM) 
MUF(I)=DARSIN(l.DO/MF(I)) 

BF(I)=0. 

3  CONTINUE 
DY=DY0 

DO  4  KC=1,KCHAR0 
CSIGNN(KC)=CSIGNF(KC) 

4  CONTINUE 

DO  5  1=1, KN 
BNCI)=0. 

5  CONTINUE 
RETURN 
END 


PRINT 


COMMON 
COMMON 
COMMON 


SUBKUUIlNb  PK1NI 

SUBROUTINE  NUMBER   3 

IMPLICIT  REAL*8CA-H,L-Z,$) 
REAL*4  XPL,YPL 

/PLUME/XPL ( 1 002, 10),YPL( 1002) 
/IPLUME/KPL,ITYPL(10) 

/VECS/XF(101),RMF(101),RPF(101),MF(101),MUF(101), 
TETAF(101),BF(101), 
*  XNC101),RMN(101),RPN(101),MN(101),MUN(101), 

>  TETAN(101),BN(101),XTEMP(101) 

COMMON/THICKY/XTHC1 002), TH( 1002) 
REAL*4  YXI,XI,XIPM,XIGRP,XIAPP,XIF 

/THICKX/YXI(20),XI(101,20),XIPM(101,20),XIGRP(101,20) 

,XIAPP(101,20),XIF(101,20) 
/GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14, 

G16,G17,G18,G19,G20 
/PAR/PAI,PAI2,DEG,XC,YC,MEXIT,MFIN,YMAX,DY0,DY,DYNEXT, 
STAB, DELTA, PSI1 , PSIF, ZETA1 , SIGMA, FRACG, EPSIL , NUO, 
TETSYM, TETLIM, DDY, DYMAX 
/STAG/RHO0,N0,P0,T0,A0,MDOTl 

/CHARAC/XCHARF(92),YCHARF(92),XCHARN(92),YCHARN(92), 
RMCARF(92),RPCARF(92),RMCARN(92),RPCARN(92), 
TCHARF(92),TCHARN(92),MUCARF(92),MUCARN(92), 
CSIGHN(92),CSIGNF(92),MCHARN(92),MCHARF(92), 
MCHARK92) 
/IPAR/JMAX,KF0,ITER0,KF,KN,IM,IP,J, 

KF2,IDEL,JDEL,JYXI,JXI,ILEAD,ILEADF,KCLEAD 
COMMON  /ROW/YF,YN,DXF,DXN 


COMMON 

'common 


COMMON 


COMMON 
COMMON 


COMMON 


10 


11 


SUM=0. 
KF1=KF-1 
DO  10  1=1, KF1 
DX=XF(I+1)-XF(I) 
GOREM=l .D0+G1*MF(I)**2 
GOREM1=1.DO+G1*MF(I+1)**2 

RATEM=RHO0*A0*MF(I   )*DSIN(TETAF( I   ) )/GOREM  X*(G6+0.5D0) 
RATEP=RHO0*A0*MF(I+l)*DSIN(TETAF(I+l))/GOREMl**(G6+0.5D0) 
SUM=SUM+DX*(RATEM+RATEP)/2.D0 
CONTINUE 
YYF=2.D0*PAI*YF 
IF(DELTA.EQ.O.)  YYF=1.D0 
MDOTFR=YYF*SUM/MDOTl 

PRINT  11,  J,KCLEAD,KF,ILEAD,YF,DY,XF(KF),MF(KF),MDOTFR 
FORMATUX, f J,KCLEAD,KF,ILEAD,YF,DY,XF(KF),MBOUND,MDOTR=f, 
1  4I5,5D12.4) 

PRINT  FLOW  FIELD  AT  Y=YF 


JET0289 
JET0290 
JET0291 
JET0292 
JET0293 
JET0294 
JET0295 
JET0296 
JET0297 
JET0298 
JET0299 
JET0300 
JET0301 
JET0302 
JET0303 
JET0304 
JET0305 
JET0306 
JET0307 
JET0308 
JET0309 
JET0310 
JET0311 
JET0312 

JET0313 

JET0314 
JET0315 
JET0316 
JET0317 
JET0318 
JET0319 
JET0320 
JET0321 
JET0322 
JET0323 
JET0324 
JET0325 
JET0326 

G15,JET0327 
JET0328 
JET0329 
JET0330 
JET0331 
JET0332 
JET0333 
JET0334 
JET0335 
JET0336 
JET0337 
JET0338 
JET0339 
JET0340 
JET0341 
JET0342 
JET0343 
JET0344 
JET0345 
JET0346 
JET0347 
JET0348 
JET0349 
JET0350 
JET0351 
JET0352 
JET0353 
JET0354 
JET0355 
JET0356 
JET0357 
JET0358 
JET0359 
JET0360 


24 


: : 


JETPR 


FORTRAN   Al 


IF(J.EQ.JMAX)  JXI=MINOCJXI,JYXI) 

IFU.EQ.l  .OR.  J.EQ.JMAX)  GO  TO  121 

IF(JXI.GT.JYXI)  GO  TO  120 

IFCYXIUXI).GT.YF+0.5D0*DY)  GO  TO  120 
121   CONTINUE 

YXI(JXI)=YF 

CALL  OPACX 
J  COMPUTE  MACH  NUMBER  FOR  CYLINDRICAL  EXPANSION  MCYL . 

F=(YF/YC)*(G7*(1.D0+G1*MEXIT**2))**G2/MEXIT 

CALL  AREAF(F,MCYL) 

PRINT  22,JXI,KCLEAD,ILEAD,KF,MCYL,YF 
22    F0RMAT(/1X,'PRINTING  NUMBER   JXI, KCLEAD, ILEAD, KF= • , 414, 
1        5X,  ,MCYL,YF=S2D14.5/) 

PRINT  1 


FORMATC/1X, 


KC 


,'      XFCI)   ', 

■   TETAF(I)   ', 

•      MF(I)   », 

•     MAB      «, 

■      MABI    • , 

•     MOBI     •, 

•    XI(I)    •, 

'    XIGRP(I)  ', 

'    XIAPP(I)  », 

'    XIPM(I)   V) 

IDEL1=1 

23 


21 
20 


1 
2 
3 

4 

IDEL1=IDEL 

IFU.EQ.l.  OR.  J.EQ.JMAX) 

DO  20  I=1,KF,IDEL1 

KC=KCLEAD+(I-ILEAD) 

IFCKC.LT. KCLEAD)  KC=0 

M0B=1.D10 

M0BI=1.D10 

MAB=1.D10 

MABI=1.D10 

MPM=MF(I) 

IF(KC.EQ.O)  GO  TO  23 

MOB=MCHARI(KC) 

IF(J.EQ.l)  GO  TO  23 

PSIPM=PAI2-DATAN((XF(I)-XC)/(YF-YO) 

ZETA=PSI1+ZETA1-PSIPM 

MPM=DSQRT( (G5*DTAN(ZETA/G5) )**2+l . DO) 

CALL  MATCHCI, MOB, MAB, MOBI, MABI) 

CONTINUE 

PRINT  21,I,KC,XF(I),TETAF(I)*DEG,MF(I),MAB,MABI,M0BI, 
1         XKI, JXI),XIGRP(I,JXI),XIAPP(I,JXI),XIPM(I,JXI) 

FORMAT(1X,2I4,10D12.4) 

CONTINUE 

120 
TO  120 


TO 
GO 


'RADIAL  MOLECULAR  THICKNESS   J 
(JJ,XTHUJ),THUJ),JJ  =  1,JMAX) 
D11.4,D10.3)) 


XTH(J),TH(J)  =  V) 


IF(J.EQ.l)  GO 

IFU.EQ.JMAX) 

JXI=JXI+1 

CONTINUE 

IF(J.LT.JMAX)  GO  TO  200 

PRINT  101 

FORMATCl') 

PRINT  102 

FORMATC1X, 

PRINT  202 

F0RMAT(/5(I5 

PRINT  101 

PRINT  103,(IPL,ITYPL(IPL),IPL=1,KPL) 

FORMATC1X, 'PLUME  TYPES   IPL , ITYPL ( IPL )=' , 
1  2(/lX,5(5X,2I4))) 

PRINT  104 

FORMATdX,  'PLUME  POINTS   J  ,  YPL  (  J  ) ,  XPL  (  J  ,  1 )  ,  XPL  (  J  ,  2) 

JDEL1=1 

DO  203  JJ=1, JMAX, JDEL1 

PRINT  204, JJ,YPL(JJ),(XPLUJ,IPL),IPL=1,KPL) 

FORMAT(1X,I5,2X,E12.4,10E11.3) 

CONTINUE 

C   WRITE  ON  TAPE7  FOR  SUBSEQUENT  PLOTTING. 
C   NO  MORE  THAN  80  CHARACTERS  PER  LINE  ON  TAPE7 . 

WRITE(7,205)  JMAX,KPL 

FORMATC8I10/8I10) 

WRITE(7,205)  ( ITYPL ( IPL ), IPL =1 , KPL ) 

DO  210  JJ=1,JMAX 

WRITE(7,211)  YPL(JJ),(XPL(JJ,IPL),IPL=1,KPL) 

F0RMAT(6E13.6/2X,6E13.6/2X,6E13.6/2X,6E13.6) 

CONTINUE 


120 


101 
102 
202 


103 


104 


204 
203 


205 


211 
210 


=  '/) 


JET0361 
JET0362 
JET0363 
JET0364 
JET0365 
JET0366 
JET0367 
JET0368 
JET0369 
JET0370 
JET0371 
JET0372 
JET0373 
JET0374 
JET0375 
JET0376 
JET0377 
JET0378 
JET0379 
JET0380 
JET0381 
JET0382 
JET0383 
JET0384 
JET0385 
JET0386 
JET0387 
JET0388 
JET0389 
JET0390 
JET0391 
JET0392 
JET0393 
JET0394 
JET0395 
JET0396 
JET0397 
JET0398 
JET0399 
JET0400 
JET0401 
JET0402 
JET0403 
JET0404 
JET0405 
JET0406 
JET0407 
JET0408 
JET0409 
JET0410 
JET0411 
JET0412 
JET0413 
JET0414 
JET0415 
JET0416 
JET0417 
JET0418 
JET0419 
JET0420 
JET0421 
JET0422 
JET0423 
JET0424 
JET0425 
JET0426 
JET0427 
JET0428 
JET0429 
JET0430 
JET0431 
JET0432 
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226 

221 
227 


JXI0,KF0=B,2I8) 


WRITE  LATERAL  (X)  OPACITIES 
JXIO=JXI 

WRITEC7,205)    JXIO,KFO 
PRINT   226,    JXICKFO 

F0RMAT(///1X,' LATERAL  (X)  OPACITIES 
DO  220  JXI=1,JXI0 
WRITE(7,221)  JXI,YXI(JXI) 
F0RMAT(I10,E13.6) 
PRINT  227,  JXI,YXI(JXI) 
F0RMAT(//1X,'JXI,YXI(JXI)=«,I8,E15.6/) 
DO  225  1=1, KFO 

WRITE(7,211)  XIF(I,JXI),XI(I,JXI),XIPM(I,JXI),XIGRP(I,JXI), 
XIAPP(I,JXI) 


PRINT  211,  XIF(I,JXI),XI(I,JXI),XIPM(I,JXI),XIGRP(I,JXI), 
1  XIAPP(I,JXI) 

225  CONTINUE 
220  CONTINUE 
200   CONTINUE 

RETURN 

END 


FIN 


JET0433 
JET0434 
JET0435 
JET0436 
JET0437 
JET0438 
JET0439 
JET0440 
JET0441 
JET0442 
JET0443 
JET0444 
JET0445 
JET0446 
JET0447 
JET0448 
JET0449 
JET0450 
JET0451 
JETQ^52 


SUBROUTINE  FINdFIN) 

C   SUBROUTINE  NUMBER   4 

C   STOP  WHEN  ERROR  IS  DETECTED. 

IMPLICIT  REAL*8(A-H,L-Z,$) 

PRINT  1,IFIN 
1     F0RMAT(/1X,»FIN  CODE  IFIN=',I6/) 
C   INDUCE  ERROR  IN  ORDER  TO  GENERATE  TRACING  OF  CALLING  SUBROUTINES. 

X=-1.D0 

Y=X+DSQRTCX) 

IF(IFIN.LE.O)  GO  TO  100 

STOP 

100    g^URN HAJ?CH 

SUBROUTINE  MARCH 
C   SUBROUTINE  NUMBER   5 

IMPLICIT  REAL*8CA-H,L-Z,$) 

COMMON  /VECS/XF(101),RMF(101),RPF(101),MF(101),MUF(101), 

1  TETAF(101),BF(101), 

2  XN(101),RMN(101),RPN(101),MN(101),MUN(101), 

3  TETAN(101),BN(101),XTEMP(101) 
COMMON/THICKY/XTH(10  02),TH(10  02) 

COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15 
1  G16,G17,G18,G19,G20 

COMMON  /PAR/PAI,PAI2,DEG,XC,YC,MEXIT,MFIN,YMAX,DY0,DY,DYNEXT, 

1  STAB, DELTA, PSI1 , PSIF, ZETA1 , SIGMA, FRACG, EPSIL , NUO , 

2  TETSYM,TETLIM,DDY,DYMAX 
COMMON  /STAG/RHO0,N0,P0,T0,A0,MDOTl 
COMMON  /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP,J, 

1  KF2,IDEL, JDEL, JYXI , JXI , ILEAD, ILEADF, KCLEAD 

COMMON  /ROW/YF,YN,DXF,DXN 
COMMON  /CHARAC/XCHARF(92),YCHARF(92),XCHARN(92),YCHARN(92), 

1  RMCARF(92),RPCARF(92),RMCARN(92),RPCARN(92), 

2  TCHARF(92),TCHARN(92),MUCARF(92),MUCARN(92), 

3  CSIGNN(92),CSIGNF(92),MCHARN(92),MCHARF(92), 

4  MCHARK92) 

COMMON  /ICHARA/KCHARP, KCHARM, KCHARO 


ADVANCE  FLOW  FIELD  FROM  YF  TO  YN 

IM  =  KF 

IP  =  KF 

YN=YF+DY 

KN=KF0 
SEMI-INVERSE  INTEGRATION  FOR  FAN  POINTS. 

CALL  SEMINV 
NEW  GRID  POINTS  (JUST  INVERSE  MARCHING). 

CALL  GRIDN 
LOAD  FLOW  VARIABLES  FROM  SEMI-INVERSE  INTEGRATION  INTO  VECTORS 

CALL  LOADC 
CHARACTERISTIC  SCHEME  INTEGRATION  FOR  INNER  POINTS  (INVERSE  MARCH) 

CALL  INVMAR 

RETURN 

END 


JET0453 
JET0454 
JET0455 
JET0456 
JET0457 
JET0458 
JET0459 
JET0460 
JET0461 
JET0462 
JET0463 
JET0464 
JET0465 
JET0466 
JET0467 
JET0468 
JET0469 
JET0470 
JET0471 
JET0472 
JET0473 
,JET0474 
JET0475 
JET0476 
JET0477 
JET0478 
JET0479 
JET0480 
JET0481 
JET0482 
JET0483 
JET0484 
JET0485 
JET0486 
JET0487 
JET0488 
JET0489 
JET0490 
JET0491 
JET0492 
JET0493 
JET0494 
JET0495 
JET0496 
JET0497 
JET0498 
JET0499 
JET0500 
JET0501 
JET0502 
JET0503 
JET0504 
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SUBROUTINE  INVMAR 
SUBROUTINE  NUMBER   6 

IMPLICIT  REAL*8(A-H,L-Z,$) 

COMMON  /VECS/XF(101),RMF(101),RPF(101),MF(101),MUFC101), 

1  TETAFC101),BFC101), 

2  XNC101),RMNU01),RPNC101),MN(101),MUN(101), 

3  TETANU01),BNC101),XTEMP<101) 
COMMON/THICKY/XTH( 1002), THC 1002) 

COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15 
1  G16,G17,G18,G19,G20 

COMMON  /PAR/PAI,PAI2,DEG,XC,YC,MEXIT,MFIN,YMAX,DY0,DY,DYNEXT, 

1  STAB, DELTA, PSI1, PSIF,ZETA1, SIGMA, FRACG, EPSIL, NUO, 

2  TETSYM,TETLIM,DDY,DYMAX 
COMMON  /STAG/RHOO,NO,PO,TO,AO,MDOT1 
COMMON  /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP,J, 

1  KF2,IDEL,JDEL,JYXI,JXI,ILEAD,ILEADF,KCLEAD 

COMMON  /ROW/YF,YN,DXF,DXN 

INTEGRATION  WITH  INVERSE  CHARACTERISTICS  FOR  NEW  P0INTCX4, Y4) . 

OLD  POINTS  ARE  (XI, YD , (X2, Y2) . 

XI  IS  OBTAINED  BY  INVERSE  C-  FROM  X4 

X2  IS  OBTAINED  BY  INVERSE  C+  FROM  X4 

NOTE  THAT  XI  MAY  BE  NEGATIVE  (E.  G.  WHEN  X4=0). 

KN1=ILEAD-1 

IF(KNl.LE.O)  CALL  FINC601) 

DO  1000  1=1, KN1 

14  =  1 

X4=XNCI) 

Y4=YN 

IF4=(IM+IP)/2 

CALL  INTERP(0,IF4,KF,X4,XF,RM4,RMF,RP4,RPF) 

CALL  RFUNC(RM4,RP4,M4,MU4,TETA4) 

M14=M4 

MU14=MU4 

TETA14=TETA4 

M24=M4 

MU24=MU4 

TETA24=TETA4 

Y1=YF 

Y2  =  YF 

Y14=(Y1+Y4)/2.D0 

Y24=(Y2+Y4)/2.D0 

X1=1.D10 

X2=1.D10 

RM4=1.D10 

RP4=1.D10 

ITER=0 

GO  TO  2 

CORRECTOR 

ITER=ITER+1 
AVERAGED  PROPERTIES  ON  C-C14) ,C+(24)  CHARACTERISTICS. 

RM14=(RM1+RM4)/2.D0 

RP14=(RP1+RP4)/2.D0 

RM24=(RM2+RM4)/2.D0 

RP24=(RP2+RP4)/2.D0 
M14,MU14,TETA14,  M24, MU24 . TETA24  AVERAGED  ON  C-,C+  CHARACTERISTICS. 

CALL  RFUNC(RM14,RP14,M14,MU14,TETA14) 

CALL  RFUNC(RM24,RP24,M24,MU24,TETA24) 

CONTINUE 
NEW  XI, X2 

X10=X1 

X20=X2 

X1=X4-DY/DTAN(TETA14-MU14) 

X2=X4-DY/DTAN(TETA24+MU24) 

IFCX2.LT.0.)  CALL  FIN(670) 

D14=DSQRT(  (X1-X4)**2+DY**2) 

D24=DSQRT((X2-X4)**2+DY**2) 
INTERPOLATE  OLD  DISTRIBUTION  FOR  RM1,RP1,  RM2,RP2  AT  XI, X2. 

CALL  INTERP(0,IM,KF,X1,XF,RM1,RMF,RP1,RPF) 

CALL  INTERP(0,IP,KF,X2,XF,RM2,RMF,RP2,RPF) 


JET0505 
JET0506 
JET0507 
JET0508 
JET0509 
JET0510 
JET0511 
JET0512 
,JET0513 
JET0514 
JET0515 
JET0516 
JET0517 
JET0518 
JET0519 
JET0520 
JET0521 
JET0522 
JET0523 
JET0524 
JET0525 
JET0526 
JET0527 
JET0528 
JET0529 
JET0530 
JET0531 
JET0532 
JET0533 
JET0534 
JET0535 
JET0536 
JET0537 
JET0538 
JET0539 
JET0540 
JET0541 
JET0542 
JET0543 
JET0544 
JET0545 
JET0546 
JET0547 
JET0548 
JET0549 
JET0550 
JET0551 
JET0552 
JET0553 
JET0554 
JET0555 
JET0556 
JET0557 
JET0558 
JET0559 
JET0560 
JET0561 
JET0562 
JET0563 
JET0564 
JET0565 
JET0566 
JET0567 
JET0568 
JET0569 
JET0570 
JET0571 
JET0572 
JET0573 
JET0574 
JET0575 
JET0576 
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C   NO  NEED  FOR  RE-AVERAGING  SINCE  IT  INTRODUCES  ONLY  HIGHER  ORDER  JET0577 

C   CHANGES  INTO  THE  ITERATION  SCHEME.  JET0578 

C   INTEGRATE  THE  CHARACTERISTIC  EQUATIONS  FOR  RM4,RP4  AT  X4,Y4.  JET0579 

RM40=RM4  JET0580 

RP40=RP4  JET0581 

RM4=RM1+DELTA*DSIN(TETA14)*D14/(M14*Y14)  JET0582 

RP4=RP2+DELTA*DSINCTETA24)*D24/CM24*Y24)  JET0583 

C   CONVERGENCE  TEST  JET0584 

EPS=(DABS(X1-X10)+DABS(X2-X20))/DY+DABSCRM4-RM40)+DABS(RP4-RP40)   JET0585 

IF(ITER.GT.ITERO)  GO  TO  10  JET0586 

IF(EPS.GT.EPSIL)  GO  TO  1  JET0587 

RMN(I)=RM4  JET0588 

RPN(I)=RP4  JET0589 

CALL  RFUNC(RM4,RP4,MN(I),MUN(I),TETAN(I))  JET0590 

1000  CONTINUE  JET0591 

RETURN  JET0592 

10  CONTINUE  JET0593 
PRINT  11,I4,KN,IF4,IM,IP,KF,ITER,ITER0,EPS,EPSIL,X1,X2,X4,M14,M24  JET0594 

11  F0RMAT(1X,«SUBR.  INVMAR.  14, KN, IF4, IM, IP, KF, ITER, ITER0= ■ ,815/  JET0595 
1        1X,'EPS,EPSIL,X1,X2,X4,M14,M24=,,7D14.6/)  JET0596 

CALL  FINC611)  JET0597 

RETURN                                          CZTM/a/I/  JET0598 

END g&MJN  V JET0599 

SUBROUTINE  SEMINV JET0600 

C   SUBROUTINE  NUMBER   7  JET0601 

IMPLICIT  REAL*8(A-H,L-Z,$)  JET0602 

COMMON  /VECS/XF(101),RMF(101),RPFC101),MF(101),MUF(101),  JET06  03 

1  TETAF(101),BF(101),  JET0604 

2  XN(101),RMN(101),RPN(101),MN(101),MUN(101),  JET0605 

3  TETAN(101),BN(101),XTEMP(101)  JET0606 
COMMON/THICKY/XTHC1 002 ),TH( 1002)  JET06  07 
COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15,JET06  08 

1              G16,G17,G18,G19,G20  JET0609 

COMMON  /PAR/PAI,PAI2,DEG,XC,YC,MEXIT,MFIN,YMAX,DY0,DY,DYNEXT,  JET0610 

1  STAB, DELTA, PSI1 , PSIF,ZETA1 , SIGMA, FRACG, EPSIL , NUO,  JET0611 

2  TETSYM,TETLIM,DDY,DYMAX  JET0612 
COMMON  /STAG/RHO0,N0,P0,T0,A0,MDOTl  JET0613 
COMMON  /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP,J,  JET0614 

1              KF2,IDEL,JDEL,JYXI,JXI,ILEAD,ILEADF,KCLEAD  JET0615 

COMMON  /ROW/YF,YN,DXF,DXN  JET0616 

COMMON  /CHARAC/XCHARF(92),YCHARF(92),XCHARN(92),YCHARN(92),  JET0617 

1  RMCARF(92),RPCARF(92),RMCARN(92),RPCARN(92),  JET0618 

2  TCHARFC92),TCHARN(92),MUCARF(92),MUCARN(92),  JET0619 

3  CSIGNN(92),CSIGNF(92),MCHARN(92),MCHARF(92),  JET0620 

4  MCHARK92)  JET0621 
COMMON  /ICHARA/KCHARP,KCHARM,KCHARO  JET0622 

C  JET0623 

C   COMPUTE  NEW  POINT  (X4,Y4),  BY  PASSING  A  C+  CHARACTERISTIC  JET0624 

C   THROUGH  OLD  POINT  (X2,Y2).   BOTH  POINTS  ARE  ON  CHARACTERISTIC  LINE    JET0625 

C   NUMBER  KC.  JET0626 

IM=1  JET0627 

DO  100  KC=1,KCHAR0  JET0628 

IF(CSIGNNCKC) . EQ . 0 . )  GO  TO  100  JET0629 

C  JET0630 

C   PREDICTOR  JET0631 

C  JET0632 

Y1=YF  JET0633 

Y2=YF  JET0634 

Y4=YN  JET0635 

Y14=(Y1+Y4)/2.D0  JET0636 

Y24=(Y2+Y4)/2.D0  JET0637 

X2=XCHARF(KC)  JET0638 

RM2=RMCARF(KC)  JET0639 

RP2=RPCARF(KC)  JET0640 

M2=MCHARF(KC)  JET0641 

MU2=MUCARF(KC)  JET0642 

TETA2=TCHARF(KC)  JET0643 

M14=M2  JET0644 

MU14=MU2  JET0645 

TETA14=TETA2  JET0646 

M24=M2  JET0647 

MU24=MU2  JET0648 
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TETA24=TETA2  JET0649 

X4=1.D10  JET0650 

X1=1.D10  JET0651 

RM4=1.D10  JET0652 

RP<+  =  1.D10  JET0653 

ITER=0  JET0654 

GO   TO    2  JET0655 

C  JET0656 

C      CORRECTOR  JET0657 

C  JET0658 

1  ITER=ITER+1                                       -  .  JET0659 
C   AVERAGED  PROPERTIES  ON  C-C14) ,C+C24)  CHARACTERISTICS.  JET0660 

RM14=(RM1+RM4)/2.D0  JET0661 

RP14=(RP1+RP4)/2.D0  JET0662 

RM24=(RM2+RM4)/2.D0  JET0663 

RP24=CRP2+RP4)/2.D0  JET0664 
C   M14,MU14,fETA14,  M24,MU24,TETA24  AVERAGED  ON  C-,C+  CHARACTERISTICS.   JET0665 

CALL  RFUNC(RM14,RP14,M14,MU14,TETA14)  JET0666 

CALL  RFUNC(RM24,RP24,M24,MU24,TETA24)  JET0667 

2  CONTINUE  JET0668 
C   NEW  X4,X1  JET0669 

X40=X4  JET0670 

X10=X1  JET0671 

X4=X2+DY/DTAN(TETA24+CSIGNF(KC)*MU24)  JET0672 

X1=X4-DY/DTAN(TETA14-CSIGNF(KC)*MU14)  JET067  3 

D14=DSQRT((X1-X4)**2+DY**2)  JET0674 

D24=DSQRT((X2-X4)**2+DY**2)  JET0675 

C   INTERPOLATE  OLD  DISTRIBUTION  FOR  RM1,RP1,  AT  XI.  JET0676 

CALL  INTERPC0,IM,KF,X1,XF,RM1,RMF,RP1,RPF)  JET0677 

IF(J.GT.l)  GO  TO  22  JET0678 

IF(CSIGNF(KC).LT.O.)  GO  TO  22  JET0679 

RP1=RP2  JET0680 

22    CONTINUE  JET0681 

C   NO  NEED  FOR  RE-AVERAGING  SINCE  IT  INTRODUCES  ONLY  HIGHER  ORDER  JET0682 

C   CHANGES  INTO  THE  ITERATION  SCHEME.  JET0683 

C   INTEGRATE  THE  CHARACTERISTIC  EQUATIONS  FOR  RM4,RP4  AT  X4,Y4.  JET0684 

RM40=RM4  JET0685 

RP40=RP4  JET0686 

IF(CSIGNFCKC) .LT.O. )  GO  TO  21  JET0687 

RM4=RM1+DELTA*DSIN(TETA14)*D14/(M14*Y14)  JET0688 

RP4  =  RP2+DELTA*DSIN(TETA24;>*D24/(M24*Y24)  JET0689 

GO  TO  20  JET0690 

21    CONTINUE  JET0691 

RM4=RM2+DELTA*DSIN(TETA24)*D24/(M24*Y24)  JET06  92 

RP4=RP1+DELTA*DSIN(TETA14)*D14/(M14*Y14)  JET06  93 

20    CONTINUE  JET0694 

C   CONVERGENCE  TEST  JET0695 

EPS=(DABS(X4-X40)+DABS(X1-X10)VDY+DABS(RM4-RM40)+DABS(RP4-RP40)   JET06  96 

IF(ITER.GT.ITERO)  GO  TO  10  JET0697 

IF(EPS.GT.EPSIL)  GO  TO  1  JET0698 

CSIGNN(KC)=CSIGNF(KC)  JET0699 

IFCX4.GT.0.)  GO  TO  30  JET0700 

RMSAVE=RM4  JET0701 

RM4=RP4+TETSYM  JET0702 

RP4=RM4-TETSYM  JET0703 

CSIGNN(KC)=-1.D0  JET0704 

30    CONTINUE  JET0705 

RMCARN(KC)=RM4  JET0706 

RPCARN(KC)=RP4  JET0707 

CALL  RFUNC(RM4,RP4,M4,MU4,TETA4)  JET0708 

TCHARN(KC)=TETA4  JET0709 

XCHARN(KC)=DABS(X4)  JET0710 

YCHARN(KC)=Y4  JET0711 

MUCARN(KC)=MU4  JET0712 

MCHARN(KC)=M4  JET0713 

100   CONTINUE  JET0714 

RETURN  JET0715 

10  CONTINUE  JET0716 
PRINT  11,KC,KCHAR0,IM,KF,ITER,ITER0,EPS,EPSIL,X1,X2,X4,M14,M24     JET0717 

11  FORMATUX, 'SUBR.  SEMINV.  KC, KCHARO , IM, KF, ITER, ITER0= ', 615/  JET0718 
1        IX, »EPS,EPSIL,X1,X2,X4,M14,M24=',7D14.6/)  JET0719 

CALL  FINC711)  JET0720 
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RETURN                                            QCriiiJr  JET0721 

.   END . Kri/ft^- ,lFTfl7?? 

SUBROUTINE  RFUNCCRM, RP,M,MU, TETA)  JET0723 

C   SUBROUTINE  NUMBER   8  JET0724 

IMPLICIT  REAL*8(A-H,L-Z,$)  JET0725 

COMMON  /VECS/XF(101),RMF(101),RPF(101),MF(101),MUF(101),  JET0726 

1  TETAF(101),BF(101),  JET0727 

2  XN(101),RMN(101),RPN(101),MNC101),MUNC101),  JET0728 

3  TETAN(101),BNC101),XTEMPC101)  JET0729 
COMMON/THICKY/XTHC 1002), TH( 1002)  JET0730 
COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15,JET0731 

1             G16,G17,G18,G19,G20  JET0732 

COMMON  /PAR/PAI,PAI2,DEG,XC,YC,MEXIT,MFIN,YMAX,DY0,DY,DYNEXT,  JET0733 

1  STAB, DELTA, PSI1 , PSIF,ZETA1 , SIGMA, FRACG, EPSIL, NUO,  JET0734 

2  TETSYM,TETLIM,DDY,DYMAX  JET0735 
COMMON  /STAG/RHO0,N0,P0,T0,A0,MDOTl  JET0736 
COMMON  /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP,J,  JET0737 

1             KF2,IDEL,JDEL,JYXI,JXI,ILEAD,ILEADF,KCLEAD  JET0738 

COMMON  /ROW/YF,YN,DXF,DXN  JET0739 

C  JET0740 

C   COMPUTE  M,MU,TETA  AT  A  POINT,  AS  FUNCTION  OF  RIEMANN  INVAR.  RM,RP.    JET0741 

TETA=(RM-RP)/2.D0+TETLIM  JET0742 

NU   =(RM+RP)/2.D0  JET0743 

C   NU=NU0-(G5*ARCTAN(G5*Q)-ARCTAN(Q)),  WHERE  Q=(M**2-1 )*x(-l/2)  JET0744 

C   FIND  Q(NU),  AND  HENCE  M(NU),  THROUGH  NEWTON  RAPHSON  ITERATIONS.  JET0745 

Q  =  -(NU-NU0)/(G4-1.D0)  JET0746 

IFCQ.LE.O.)  CALL  FINC801)  JET0747 

ITER=0  JET0748 

1     ITER=ITER+1  JET0749 

QF=Q  JET0750 

DNUDT=-(G4-1.D0)/C(1.D0+G4*QXX2)*(1.D0+Q**2))  JET0751 

DNU=NU-(NU0-(G5*DATAN(G5*Q)-DATAN(Q)))  JET0752 

Q=Q+DNU/DNUDT  JET0753 

IFCQ.LE.O.)  CALL  FINC811)  JET0754 

EPS=DABS(Q-QF)/Q  JET0755 

IF(ITER.GT.ITERO)  GO  TO  10  JET0756 

IFCEPS.GT.EPSILXl.D-3)  GO  TO  1  JET0757 

M=DSQRT(1.D0+1.D0/Q**2)  JET0758 

MU=DARSIN(1.D0/M)  JET0759 

RETURN  JET0760 

10    CONTINUE  JET0761 

CALL  FINC810)  JET0762 

RETURN                                            IM-r£QP  JET0763 

END     (fV      ' JET0764 

SUBROUTINE  INTERPUNEW,  I ,  KGRID,  X,  XVEC,  RM,  RMVEC,  RP,  RPVEC)  JET0765 

C   SUBROUTINE  NUMBER   9  JET0766 

IMPLICIT  REAL*8(A-H,L-Z,$)  JET0767 

COMMON  /VECS/XF(101),RMF(101),RPF(101),MF(101),MUF(101),  JET0768 

1  TETAF(101),BF(101),  JET0769 

2  XN(101),RMN(101),RPN(101),MN(101),MUN(101),  JET0770 

3  TETAN(101),BN(101),XTEMP(101)  JET0771 
COMMON/THICKY/XTH(10  02),TH(1002)  JET0772 
COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15,JET077  3 

1              G16,G17,G18,G19,G20  JET0774 

COMMON  /PAR/PAI,PAI2,DEG,XC,YC,MEXIT,MFIN,YMAX,DY0,DY,DYNEXT,  JET0775 

1  STAB, DELTA, PSI1 , PSIF, ZETA1 , SIGMA, FRACG, EPSIL , NUO ,  JET0776 

2  TETSYM,TETLIM,DDY,DYMAX  JET0777 
COMMON  /STAG/RHOO,NO,PO,TO,AO,MDOT1  JET0778 
COMMON  /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP,J,  JET0779 

1              KF2,IDEL,JDEL,JYXI,JXI,ILEAD,ILEADF,KCLEAD  JET0780 

COMMON  /ROW/YF,YN,DXF,DXN  JET0781 

DIMENSION  XVEC(1),RMVEC(1),RPVEC(1)  JET0782 

C  JET0783 

C   FIND  I  SUCH  THAT  XVEC( I ) . LE . X . AND . XVECC 1+1 ) . GE.X  JET0784 

C   FIND  RM,RP  BY  LINEAR  INTERPOLATION.  JET0785 

C   NOTE  THAT  X  MAY  BE  NEGATIVE.  JET0786 

IF(DABS(X) .LE.XVEC(KGRID))  GO  TO  901  JET0787 

PRINT  900,  X,  KGRID, XVECCKGRID)  JET0788 

FORMAT(/1X,D15.7,I10,4X,D15.7/)  JET0789 

CALL  FIN(900)  JET0790 

901   CONTINUE  JET0791 

KG2=2*KGRID  JET0792 
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IO=MINO(I,KGRID-2) 
ICOUNT=0 

I  1  =  10 
SIGN1=1.D0 
IF(I.GE.l)  GO  TO  10 
SIGN1=-1.D0 
I=-I+2 

10    CONTINUE 

IF(I.GT.KGRID)  CALL  FINC901) 

XX1=SIGN1*XVEC(I) 

11  =  1 

IF(XXl.LE.X)  GO  TO  11 

10=10-1 

ICOUNT=ICOUNT+l 

IFCIC0UNT.GT.KG2)  CALL  FINC911) 

GO  TO  1 

II  CONTINUE 
1=10+1 
SIGN2=1.D0 
IF(I.GE.l)  GO  TO  12 
SIGN2=-1.D0 
I=-I+2 

12  CONTINUE 
IF(I.GT.KGRID)  CALL  FIN(912) 
XX2=SIGN2*XVECCI) 

12=1 

IFCXX2.GE.X)  GO  TO  13 

10=10+1 

ICOUNT=ICOUNT+l 

IFCIC0UNT.GT.KG2)  CALL  FINC913) 

GO  TO  1 

13  CONTINUE 
F1=(XX2-X)/(XX2-XX1) 
F2=1.D0-F1 

IFCF1.LT.0.)  CALL  FINC991) 
IFCF2.LT.0.)  CALL  FINC992) 
RM1=RMF(I1) 
RP1=RPF(I1) 
RM2=RMF(I2) 
RP2=RPF(I2) 
IFCSIGN1.LT.0 
IFCSIGNl.LT. 0 
IF(SIGN2.LT.O 
IFCSIGN2.LT.0 
RM=F1*RM1+F2*RM2 
RP=F1*RP1+F2*RP2 
RETURN 
END 


)  RM1=RPF(I1)+TETSYM 

)  RP1=RMF(I1)-TETSYM 

)  RM2=RPF(I2)+TETSYM 

)  RP2=RMF(I2)-TETSYM 


ItJTEQX 


bUBRUUIlNh    lNlbKXUNbW,ll,VAKU,VAK,KGKlU,XU,XVbC) 

SUBROUTINE  NUMBER   10 

IMPLICIT  REAL*8(A-H,L-Z,$) 

COMMON  /VECS/XF(101),RMF(101),RPF(101),MF(101),MUF(101) 

1  TETAF(101),BF(101), 

2  XN(101),RMN(101),RPN(101),MN(101),MUN(101) 

3  TETAN(101),BN(101),XTEMP(101) 
COMMON/THICKY/XTH(1002),TH(10  02) 


JET0793 
JET0794 
JET0795 
JET0796 
JET0797 
JET0798 
JET0799 
JET0800 
JET0801 
JET0802 
JET0803 
JET08Q4 
JET0805 
JET0806 
JET0807 
JET0808 
JET0809 
JET0810 
JET0811 
JET0812 
JET0813 
JET0814 
JET0815 
JET0816 
JET0817 
JET0818 
JET0819 
JET0820 
JET0821 
JET0822 
JET0823 
JET0824 
JET0825 
JET0826 
JET0827 
JET0828 
JET0829 
JET0830 
JET0831 
JET0832 
JET0833 
JET0834 
JET0835 
JET0836 
JET0837 
JET0838 
JET0839 
JET0840 


JET0841 
JET0842 
JET0843 
JET0844 
JET0845 
JET0846 
JET0847 
JET0848 


COMMON  /GAMA/G,G1,G2,G3,G<+,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15,JET0849 

1             G16,G17,G13,G19,G20  JET0850 

COMMON  /PAR/PAI,PAI2,DEG,XC,YC,MEXIT,MFIN,YMAX,DY0,DY,DYNEXT,  JET0851 

1  STAB,DELTA,PSI1,PSIF,ZETA1,SIGMA,FRACG,EPSIL,NU0>  JET0852 

2  TETSYM,TETLIM,DDY,DYMAX  JET0853 
COMMON  /STAG/RH00,N0,P0,T0,A0,MD0T1  JET0854 
COMMON  /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP, J,  JET0855 

1             KF2,IDEL,JDEL,JYXI,JXI,ILEAD,ILEADF,KCLEAD  JET0856 

COMMON  /ROW/YF,YN,DXF,DXN  JET0857 

DIMENSION  VAR(1),XVEC(1)  JET0858 

FIND  XO  AND  II  SUCH  THAT  XVEC( II )<X0<XVEC( 11+1 ) ,  AND  XO  CORRESPONDS  JET0859 

TO  THE  LOCATION  AT  WHICH  VARO  IS  A  LINEAR  INTERPOLATION  OF  VARCI).  JET0860 

X0=1.D23  JET0861 

IFIRST=I1  JET0362 

IF(Il.GT.O)  GO  TO  10  JET0863 

IFIRST=KGRID-IABS(Il)+2  JET0864 
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10  CONTINUE 

DO  1  II=IFIRST,KGRID 

1  =  11 

IF(Il.GT.O)  GO  TO  11 

I=KGRID-II+2 

11  CONTINUE 

IFCI.LE.O)  CALL  FIN(lOOl) 

IF(I.GT.KGRID)  CALL  FINC1002) 

IFCI.EQ.l)  GO  TO  1 

IF((VAR(I)-VAR0)*(VARCI-1)-VAR0).GT.0.) 

IFCVARCI).EQ.VAR(I-l))  GO  TO  1 

F1=CVAR(I)-VAR0)/(VARCI)-VAR(I-1)) 

F2=1.D0-F1 

IFCF1.LT.0.)  CALL  FINC1011) 

IFCF2.LT. 0.)  CALL  FINU012) 

X0=F1*XVEC(I-1)+F2*XVEC(I) 

11=1-1 

GO  TO  2 

1  CONTINUE 

2  CONTINUE 
RETURN 
END 


GO  TO  1 


B1ZSAK 


— SUBKUUIlNh  BRhAK 

C   SUBROUTINE  NUMBER   11 

IMPLICIT  REALX8(A-H,L-Z,$) 

REAL  MB, MX, MY 

COMMON  /VECS/XF(101),RMF(101),RPF(101),MF(101),MUF(101), 

1  TETAF(101),BF(101), 

2  XN(101),RMN(101),RPN(101),MN(101),MUN(101), 

3  TETAN(101),BN(101),XTEMP(101) 
COMMON/THICKY/XTHU 002 ),TH( 1002) 

COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15 
1  G16,G17,G18,G19,G20 

COMMON  /PAR/PAI,PAI2,DEG,XC,YC,MEXIT,MFIN,YMAX,DY0,DY,DYNEXT, 

1  STAB, DELTA, PSI1 , PSIF, ZETA1 , SIGMA, FRACG, EPSIL, NUO, 

2  TETSYM,TETLIM,DDY,DYMAX 
COMMON  /STAG/RHO0,N0,P0,T0,A0,MDOTl 
COMMON  /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP,J, 

1  KF2,IDEL,JDEL,JYXI,JXI,ILEAD,ILEADF,KCLEAD 

COMMON  /ROW/YF,YN,DXF,DXN 


JET0865 
JET0866 
JET0867 
JET0868 
JET0869 
JET0870 
JET0871 
JET0872 
JET0873 
JET0874 
JET0875 
JET0876 
JET0877 
JET0878 
JET0879 
JET0880 
JET0881 
JET0882 
JET0883 
JET0884 
JET0885 
JET0886 


COMPUTE  THE  BREAKDOWN  PARAMETER  AT  ( 1-1/2, K-l/2) . 

YB=0.5D0*(YF+YN) 

DYY=DY 

IM  =  2 

DO  1  1=2, KN 

X1=XN(I-1) 

X2=XN(I) 

DXX=X2-X1 

IF(X2.GT.XF(KF))  GO  TO  2 

CALL  INTERP(0,IM,KF,X1,XF,RM1,RMF,RP1,RPF) 
INTERP(0,IM,KF,X2,XF,RM2,RMF,RP2,RPF) 
RFUNC(RM1,RP1,M1,MU1,TETA1) 
RFUNC(RM2,RP2,M2,MU2,TETA2) 
5D0*((MN(I)-MN(I-1))+(M2-M1))/DXX 
5D0*((MN(I)-M2)+(MN(I-1)-M1))/DYY 
25D0*(MN(I-1)+MN(I)+M1+M2) 
:Q.25D0*(TETAN(I-1)+TETAN(I)+TETA1+TETA2) 


STORE  IN  BN(I) 


CALL 

CALL 

CALL 

MX  =  0 

MY  =  0 

MB  =  0 

TETAB  = 

G0REM=MB**2*(1.D0+G1*MB**2)**(G6-1 

GRAD=MX*DCOSCTETAB)+MY*DSIN(TETAB) 

B=G20*GOREM*GRAD 

GO  TO  3 

B=1.D22 

BN(I)=B 

CONTINUE 

BN(1)=BN(2) 

RETURN 

END 


DO) 


oPAcy 


JETTJSS7 
JET0888 
JET0889 
JET0890 
JET0891 
JET0892 
JET0893 
JETQ894 
JET0895 
,JET0896 
JET0897 
JET0898 
JET0899 
JET0900 
JET0901 
JET0902 
JET0903 
JET0904 
JET0905 
JET0906 
JET0907 
JET0908 
JET0909 
JET0910 
JET0911 
JET0912 
JET0913 
JET0914 
JET0915 
JET0916 
JET0917 
JET0918 
JET0919 
JET0920 
JET0921 
JET0922 
JET0923 
JET0924 
JET0925 
JET0926 
JET0927 
JET0928 
JET0929 
JET0930 
JET0931 
JET0932 


SUBKUUIlNb  UHACY 

SUBROUTINE  NUMBER   12 

IMPLICIT  REAL*8(A-H,L-Z,$) 

COMMON  /VECS/XF(101),RMF(101),RPF(101),MF(101),MUF(101) 


Jhl U933 

JET0934 
JET0935 
JET0936 
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1 

11 


1  TETAF(101),BF(101), 

2  XN(101),RMN(101),RPN(101),MN(101),MUNC101), 

3  TETAN(101),BNC101),XTEMP(101) 
COMMON/THICKY/XTHU 002), TH( 1002) 

COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15, 
1  G16,G17,G18,G19,G20 

COMMON  /PAR/PAI,PAI2,DEG,XC,YC,MEXIT,MFIN,YMAX,DY0,DY,DYNEXT, 

1  STAB, DELTA, PSI1 , PSIF, ZETA1 , SIGMA, FRACG, EPSIL , NUO, 

2  TETSYM,TETLIM,DDY,DYMAX 
COMMON  /STAG/RHO0,N0,P0,T0,A0,MDOTl 
COMMON  /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP,J, 

1  KF2,IDEL,JDEL,JYXI,JXI,ILEAD,ILEADF,KCLEAD 

COMMON  /ROW/YF,YN,DXF,DXN 

COMPUTE  THE  MOLECULAR  THICKNESS  AT  END  POINTS  OF  EACH  ROW. 
IM=2 

XTH(J)=XF(KF) 
TH(J)=0. 

DTHO=NO*SIGMA*DY 
IF(J.EQ.l)  GO  TO  11 
J1=J-1 

DO  1  JJ=1,J1 
XXI =XTH( J J ) 

CALL  INTERP(0,IM,KF,XX1,XF,RM1,RMF,RP1,RPF) 
CALL  RFUNC(RM1,RP1,M1,MU1,TETA1) 
GOREM=1.DO+G1*M1**2 
DTH=DTH0/GOREM**G6 
TH(JJ)=TH(JJ)+DTH 
CONTINUE 
CONTINUE  _ 

*^URN  ?LUM£S 


SUBKUUIlNb  PLUMES 

SUBROUTINE  NUMBER   13 

IMPLICIT  REAL*8(A-H,L-Z,$) 

REALX4  XPL  YPL 

COMMON  /PLUME/XPLC 1002, 10), YPL (1002) 

COMMON  /IPLUME/KPL,ITYPL(10) 

DIMENSION  VPLC92) 

COMMON  /VECS/XF(101),RMF(101),RPF(101),MF(101),MUF(101) 

1  TETAF(101),BF(101), 

2  XN(101),RMN(101),RPNC101),MN(101),MUN(101) 

3  TETAN(101),BN(101),XTEMP(101) 
COMMON/THICKY/XTHC1 002 ),TH( 1002) 
REAL*4  YXI,XI,XIPM,XIGRP,XIAPP,XIF 
COMMON  /THICKX/YXI(20),XI(101,20),XIPM(101,20) 

1  ,XIAPP(101,20),XIF(101,20) 

COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15, 
1  G16,G17,G13,G19,G20 

COMMON  /PAR/PAI,PAI2,DEG,XC,YC,MEXIT,MFIN,YMAX,DY0,DY,DYNEXT, 
1  STAB, DELTA, PSI1 , PSIF, ZETA1 , SIGMA, FRACG, EPSIL, NUO, 

TETSYM,TETLIM,DDY,DYMAX 
/STAG/RHOO,NO,PO,TO,AO,MDOT1 
/IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP,J, 

KF2,IDEL, JDEL,JYXI,JXI,ILEAD,ILEADF,KCLEAD 
/ROW/YF,YN,DXF,DXN 
/CHARAC/XCHARF(92),YCHARF(92) 
RMCARF(92),RPCARF(92) 
TCHARF(92),TCHARN(92) 
CSIGNN(92),CSIGNF(92),MCHARN(92) 
MCHARK92) 
COMMON  /ICHARA/KCHARP, KCHARM, KCHARO 
COMPUTE  SPECIAL  POINTS  AT  Y=YN,  AND  STORE  THEM  AS 

(XPL(J,IPL),YPL(J)=YN). 
J  IS  THE  MARCHING  INDEX  OF  YN . 

IPL=1,2, . . .,KPL  IS  THE  "PLUME"  INDEX.   PRESENTLY   KPL.LE.5 
VPL(IPL)  IS  A  VALUE  DEFINING  THE  "PLUME"  CURVE. 

THE  TYPE  OF  CURVE.  IT  DEFINES  CURVES  AS  FOLLOWS: 
DO  NOTHING 

REAL  PLUME.  IT  IS  THE  BREAKDOWN  SURFACE,  DEFINED 
BY  A  CONSTANT  VALUE  OF  THE  BREAKDOWN  PARAMETER  B. 
SET   VPL(IPL)=B. 


JET0937 
JET0938 
JET0939 
JET0940 
JET0941 
JET0942 
JET0943 
JET0944 
JET0945 
JET0946 
JET0947 
JET0948 
JET0949 
JET0950 
JET0951 
JET0952 
JET0953 
JET0954 
JET0955 
JET0956 
JET0957 
JET0958 
JET0959 
JET0960 
JET0961 
JET0962 
JET0963 
JET0964 
JET0965 
JET0966 
JET0967 
JET0968 


JET0969 
JET0970 
JET0971 
JET0972 
JET0973 
JET0974 
JET0975 
JET0976 
JET0977 
JET0978 
JET0979 
JET0980 
JET0981 
JET0982 
JET0983 
JET0984 
JET0985 
JET0986 
JET0987 
JET0988 
JET0989 
JET0990 
JET0991 
JET0992 
JET0993 
JET0994 
JET0995 
JET0996 
JET0997 
JET0998 
JET0999 
JET1000 
JET1001 
JET1002 
JET1003 
JET1004 
JET1005 
JET1006 
JET1007 
JET1008 


XIGRP(101,20) 


COMMON 
COMMON 

COMMON 
COMMON 


XCHARN(92),YCHARN(92) 
RMCARN(92),RPCARNC92) 
MUCARF(92),MUCARN(92) 
MCHARF(92) 


C 
C 
C 
C 
C 
C 
C 
C 

c 
c 


ITYPL(IPL)  IS 

ITYPL(IPL)=0 

ITYPL(IPL)=1 
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C   ITYPL(IPL)=2   CONSTANT  MACH-NUMBER 

LINE. 

VPL(IPL)=M. 

JET1009 

C   ITYPL(IPL)=3   A  SINGLE  STREAMLINE. 

VPL(IPL) 

IS  SET  TO 

THE  EXI' 

JET1010 

C 

X-COORDINATE  OF  THAT 

STREAMLINE 

JET1011 

C   ITYPL(IPL)=4   A  SINGLE  C+  CHARACTERISTIC 

LINE 

^STARTING 

AT  THE 

:  CORNER 

.JET1012 

C 

VPL(IPL)  IS  SET  TO  THE  INDEX  KC 

:  OF  THAT  CH 

JET1013 

C 

LINE. 

JET1014 

C   ITYPL(IPL)=5   A  CONSTANT  LATERAL  (X)  OPACITY 

LINE.   VPL( 

IS  SET 

JET1015 

C 

TO  THE  VALUE  OF  THE  I 

[CONSTANT 

OPACITY. 

JET1016 

C 

JET1017 

C   DEFII                  PLC  PL) 

JET1018 

KPL=10 

JET1019 

IF(KPL.GT.IO)  CALL  FINC1301) 

JET1020 

DO  2000  IPL=1,KPL 

JET1021 

GO  TO  (2001,2002,2003,2004,2005 

,2006,20 

2008,2009, 

.2010), 

IPL 

JET1022 

2001 

ITYPL(IPL)=4 
VPL(IPL)=1 
GO  TO  2000 

JET1023 
JET1024 
JET1025 

2002 

ITYPL(IPL)=4 
VPL(IPL)=KCHARP 
GO  TO  2000 

JET1026 
JET1027 
JET1028 

2003 

ITYPL(IPL)=4 
VPL(IPL)=19 
GO  TO  2000 

JET1029 
JET1030 
JET1031 

2004 

ITYPL(IPL)=4 
VPL(IPL)=31 
GO  TO  2000 

JET1032 
JET1033 
JET1034 

2005 

ITYPL(IPL)=4 
VPL(IPL)=47 
GO  TO  2000 

JET1035 
JET1036 
JET1037 

2006 

ITYPL(IPL)=4 
VPL(IPL)=55 
GO  TO  2000 

JET1038 
JET1039 
JET1040 

2007 

ITYPL(IPL)=1 
VPL(IPL)=0.02D0 
GO  TO  2000 

JET1041 
JET1042 
JET1043 

2008 

ITYPL(IPL)=1 
VPL(IPL)=0.03D0 
GO  TO  2000 

JET1044 
JET1045 
JET1046 

2009 

ITYPL(IPL)=1 
VPL(IPL)=0.05D0 
GO  TO  2000 

JET1047 
JET1048 
JET1049 

2010 

ITYPL(IPL)=1 
VPL(IPL)=0.08D0 
GO  TO  2000 

JET1050 
JET1051 
JET1052 

2000 

CONTINUE 

JET1053 

C   COMPUT 

JET1054 

DO  1000  IPL=1,KPL 

JET1055 

ITYP=ITYPL(IPL) 

JET1056 

IF(ITYP.EQ.O)  GO  TO  1000 

JET1057 

GO  TO  (1,2,3,4,5),  ITYP 

JET1058 

1 

CONTINUE 

JET1059 

C   BREAKDOWN  SURFACE  PLUME. 

JET1060 

C   NOTE  THAT  DUE  TO  DIFFERENCE-CENTERING  OF 

GRADI 

JET1061 

C   Y- 

COORDINATE  IS  0.5*(YF+YN),  RATHER  THAN 

YN. 

IT  CAN  BE  ADJUSTED 

JET1062 

C   IN 

THE  PLOTTING  CODE. 

B0=VPL(IPL) 

XTEMP(1)=XN(1) 

DO  11  1=2, KN 

XTEMP(I)=0.5D0*(XN(I)+XN(I-1)) 

JET1063 
JET1064 
JET1065 
JET1066 
JET1067 

11 

CONTINUE 
1  =  2 

JET1068 
JET1069 

CALL  INTERX(l,I,B0,BN,KN,XBO,XTEMP) 

JET1070 

XPL(J,IPL)=XBO 

JET1071 

GO  TO  1001 

JET1072 

2 

CONTINUE 

JET1073 

c  fin: 

M=MPL. 

JET1074 

IF(J.GT.l)  GO  TO  200 

JET1075 

XPL(J,IPL)=XC 

JET1076 

GO  TO  1001 

JET1077 

200 

CONTINUE 
MPL=VPL(IPL) 
I  =  -KN 

JET1078 
JET1079 
JET1080 
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CALL  INTERX(l,I,MPL,MN,KN,XMO,XN)  JET1081 

XPL(J,IPL)=XMO  JET1082 

GO  TO  1001  JET1083 

3  CONTINUE  JET1084 
C   STREAMLINE  INTERPOLATION.  JET1085 

IF(J.GT.l)  GO  TO  300  JET1086 

XPL(J,IPL)=VPLCIPL)  JET1087 

GO  TO  1001  JET1088 

300  CONTINUE  JET1089 
XSF=XPLCJ-1,IPL)  JET1090 
ISF=2  JET1091 
ISN=2  JET1092 
CALL  INTERP(0,ISF,KF,XSF,XF,RMSF,RMF,RPSF,RPF)  JET1093 
CALL  RFUNC(RMSF,RPSF,MSF,MUSF,TETASF)  JET1094 
XSN=XSF+DY*DTAN(PAI2-TETASF)  JET1095 
ITER=1  JET1096 

301  ITER=ITER+1  JET1097 
CALL  INTERP(1,ISN,KN,XSN,XN,RMSN,RMN,RPSN,RPN)  JET1098 
CALL  RFUNC(RMSN,RPSN,MSN,MUSN,TETASN)  JET1099 
TETAAV=0.5D0*(TETASF+TETASN)  JET1100 
XSN=XSF+DY*DTAN(PAI2-TETAAV)  JET1101 
IFCITER.LT. ITER0+2)  GO  TO  301  JET1102 
XPL(J,IPL)=XSN  JET1103 
GO  TO  1001  JET1104 

4  CONTINUE  JET1105 
C   CHARACTERISTIC  LINE.  JET1106 

KC=IDINTCVPLCIPL)+l.D-5)  JET1107 

IF(J.GT.l)  GO  TO  41  JET1108 

XPL(J,IPL)=XCHARF(KC)  JET1109 

GO  TO  1001  JET1110 

41    CONTINUE  JET1111 

XPL(J,IPL)=XCHARN(KC)  JET1112 

IFCCSIGNNCKO.EQ.O.)  XPL( J, IPL)=1 . E33  JET1113 

GO  TO  1001  JET1114 

5  CONTINUE  JET1115 
C   CONSTANT  LATERAL  (X)  OPACITY  JET1116 

CALL  OPACX  JET1117 

XIC=VPL(IPL)  JET1118 

DO  51  11=2, KF  JET1119 

I1=KF-II+1  JET1120 

12=11+1  JET1121 

XI1=XI(I1,JXI)  JET1122 

XI2=XI(I2, JXI)  JET1123 

IF((XIC-XI1)*(XIC-XI2).GT.0.)  GO  TO  51  JET1124 

F2=(XI2-XIC)/(XI2-XI1)  JET1125 

F1=1.D0-F2  JET1126 

IFCF1.LT.0.)  CALL  FIN(1351)  JET1127 

IFCF2.LT.0.)  CALL  FINC1352)  JET1128 

XIFC=F2*XF(I1)+F1*XF(I2)  JET1129 

GO  TO  52  JET1130 

51  CONTINUE  JET1131 
XIFC=1.D30  JET1132 

52  CONTINUE  JET1133 
XPL(J,IPL)=XIFC  JET1134 
GO  TO  1001  JET1135 

1001  CONTINUE  JET1136 
IF(J.GT.l)  GO  TO  1002  JET1137 
YPL(J)=YC  JET1138 
GO  TO  1000  JET1139 

1002  CONTINUE  JET1140 
YPL(J)=YN  JET1141 

1000  CONTINUE  JET1142 

RETURN                                             rOllSM  JET1143 

END                                                (j-KUSW  JET1144 

SUBKUUIlNb  (JK1DN Jtlll4b 

C   SUBROUTINE  NUMBER   14  JET1146 

IMPLICIT  REAL*8(A-H,L-Z,$)  JET1147 

REAL*4  XPL,YPL  JET1148 

COMMON  /PLUME/XPL(1002,10),YPL(1002)  JET1149 

COMMON  /IPLUME/KPL,ITYPL(10)  JET1150 

COMMON  /VECS/XF(101),RMFC101),RPF(101),MF(101),MUF(101),  JET1151 

1              TETAF(101),BF(101),  JET1152 
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2  XN(101),RMNC101),RPN(101),MNC101),MUN(101),  JET1153 

3  TETAN(101),BN(101),XTEMP(101)  JET1154 
COMMON/THICKY/XTHU 002), TH( 1002)  JET1155 
COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15,JET1156 

1              G16,G17,G18,G19,G20  JET1157 

COMMON  /PAR/PAI,PAI2,DEG,XC,YC,MEXIT,MFIN,YMAX,DY0,DY,DYNEXT,  JET1158 

1  STAB, DELTA, PSI1, PSIF,ZETA1, SIGMA, FRACG, EPSIL , NUO,  JET1159 

2  TETSYM,TETLIM,DDY,DYMAX  JET1160 
COMMON  /STAG/RHO0,N0,P0,T0,A0,MDOTl  JET1161 
COMMON  /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP,J,  JET1162 

1             KF2,IDEL,JDEL,JYXI,JXI,ILEAD,ILEADF,KCLEAD  JET1163 

COMMON  /ROW/YF,YN,DXF,DXN  JET1164 

COMMON  /CHARAC/XCHARF(92),YCHARF(92),XCHARN(92),YCHARN(92),  JET1165 

1  RMCARF(92),RPCARF(92),RMCARN(92),RPCARN(92),  JET1166 

2  TCHARF(92),TCHARN(92),MUCARF(92),MUCARN(92),  JET1167 

3  CSIGNN(92),CSIGNF(92),MCHARN(92),MCHARF(92),  JET1168 

4  MCHARK92)  JET1169 
COMMON  /ICHARA/KCHARP,KCHARM,KCHARO  JET1170 

C   DIVIDE  LINE  Y=YN  INTO  KN-1  INTERVALS.  JET1171 

C   THE  X-GRID  IS  NON-UNIFORMLY  DEFINED  AS  FOLLOWS:  JET1172 

C   (1)  (XCHARN(I),YCHARN(D),  (XCHARFC I ) , YCHARFC I ) ) ,  1=1 , 2, . . . , KCHARP,   JET1173 

C       DENOTE  NEW  AND  OLD  (FORMER)  CHARACTERISTIC  (C+)  POINTS.   LET  1=1  JET1174 

C       AND  I=KCHARP  CORRESPOND  TO  THE  LEADING  AND  BOUNDARY  JET1175 

C       CHARACTERISTICS  (C+).  JET1176 

C   (2)  THE  GRID  CONSISTS  OF  TWO  SEGMENTS.   THE  SO-CALLED  FLAT  SEGMENT  JET1177 

C       IS  BETWEEN  X=0  AND  X=XLEAD=XCHARN(KCLEAD) .   THE  SECOND  IS  THE  JET1178 

C       FAN  SEGMENT.   IT  IS  FROM  XLEAD  TO  XBOUND=XCHARN(KCHARP) .  JET1179 

C   (3)  THE  FAN  SEGMENT  IS  INITIALLY  DIVIDED  INTO  FRACG*(KF0-1)  INTERVALSJET1180 

C       DEFINED  BY  THE  FAMILY  OF  C+  CHARACTERISTIC  LINES  MCHARI(l)  TO  JET1181 

C       MCHARKKCHARP).  JET1182 

C   (4)  THE  FLAT  SEGMENT  IS  DIVIDED  INTO  ( 1-FRACG)*(KF0-1 )  EQUAL  JET1183 

C       INTERVALS,  AS  LONG  AS  THEY  ARE  NOT  SMALLER  THAN  THE  AVERAGE  JET1184 

C       FAN  INTERVAL.   WHEN  THEY  ARE,  THEIR  NUMBER  IS  REDUCED,  BUT  NOT  JET1185 

C       BELOW  THREE.  JET1186 

C   (5)  KCLEAD  IS  INITIALLY  1.   IT  IS  UPDATED  SO  THAT  THE  FLAT  SEGMENT  JET1187 

C      IS  AT  LEAST  TWICE  THE  AVERAGE  FAN  INTERVAL.  JET1188 

ILEADF=ILEAD  JET1189 

KCLEAD=0  JET1190 

DO  1  KC=1, KCHARP  JET1191 

IF(CSIGNNCKC) .LT.O. )  GO  TO  1  JET1192 

KCLEAD=KC  JET1193 

KFAN=KCHARP-KCLEAD  JET1194 

XLEAD=XCHARN(KCLEAD)  JET1195 

XBOUND=XCHARN(KCHARP)  JET1196 

DX1  =  (XBOUND-XLEADVDFLOAT(KFAN)  JET1197 

IF(XLEAD/DX1.GT.2.D0)  GO  TO  11  JET1198 

1  CONTINUE  JET1199 
11    CONTINUE  JET1200 

IFCKCLEAD.EQ.      0)  CALL  FINC1401)  JET1201 

IFCKCLEAD.EQ. KCHARP)  CALL  FINU402)  JET1202 

ILEAD=IDINT(XLEAD/DXl)+2  JET1203 

IFCILEAD+KFAN.GT.KFO)  ILEAD=KFO-KFAN  JET1204 

ILEAD1=ILEAD-1  JET1205 

KN=ILEAD+KFAN  JET1206 

IF(KN.GT.KFO)  CALL  FIN(1411)  JET1207 

DX=XLEAD/DFL0AT(ILEAD1)  JET1208 

XN(1)=0.  JET1209 

DO  2  I=1,ILEAD1  JET1210 

XN(I)=XN(l)+DX*DFLOAT(I-l)  JET1211 

2  CONTINUE  JET1212 
DO  3  I=ILEAD,KN  JET1213 
XN(I)=XCHARN(KCLEAD+I-ILEAD)  JET1214 

3  CONTINUE  ^  JET1215 
RETURN                                            M^T"£P  JET1216 

END ^°  JET1217 

<jUJiR0UTINE  y5TEp Jbllk!18 

C   SUBROUTINE  NUMBER   15  JET1219 

IMPLICIT  REAL*8(A-H,L-Z,$)  JET1220 

REAL*4  XPL,YPL  JET1221 

COMMON  /PLUME/XPLU002,10),YPL(1002)  JET1222 

COMMON  /IPLUME/KPL,ITYPL(10)  JET1223 

COMMON  /VECS/XF(101),RMF(101),RPF(101),MF(101),MUF(101),  JET1224 
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1  TETAF(101),BFC101), 

2  XN(101),RMN(101),RPN(101),MN(101),MUN(101), 

3  TETAN(101),BN(101),XTEMPC101) 
COMMON/THICKY/XTHU 002), TH( 1002) 


COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15,J 

1              G16,G17,G18,G19,G20  J 

COMMON  /PAR/PAI,PAI2,DEG,XC,YC,MEXIT,MFIN,YMAX,DY0,DY, DYNEXT,  J 

1  STAB, DELTA, PSI1 , PSIF,ZETA1, SIGMA, FRACG, EPSIL, NUO,  J 

2  TETSYM,TETLIM,DDY,DYMAX  J 
COMMON  /STAG/RHO0,N0,P0,T0,A0,MDOTl  J 
COMMON  /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP,J,  J 

1              KF2,IDEL,JDEL,JYXI,JXI,ILEAD,ILEADF,KCLEAD  J 

COMMON  /ROW/YF,YN,DXF,DXN  J 

COMMON  /CHARAC/XCHARF(92),YCHARF(92),XCHARN(92),YCHARN(92),  J 

1  RMCARF(92),RPCARF(92),RMCARN(92),RPCARN(92),  J 

2  TCHARF(92),TCHARN(92),MUCARF(92),MUCARN(92),  J 

3  CSIGNN(92),CSIGNF(92),MCHARN(92),MCHARF(92),  J 

4  MCHARK92)  J 
COMMON   /ICHARA/KCHARP,KCHARM,KCHARO  J 

COMPUTE  NEXT  Y-STEP.  J 
DYNEXT  IS  DEFINED  AS  THE  MINIMAL  "TRIANGULATION"  Y-STEP  DYT,  OBTAINEDJ 

BY  FORWARD  INTERSECTION  OF  C-,C+  CHARACTERISTICS  FROM  ADJACENT  GRID  J 


POINTS  XI, X2. 

DYMIN=1.D40 

DO  1  1=3, KF 

X1=XF(I-1) 

X2=XF(I) 

DX=X2-X1 

TP1  =  DTAN(TETAF(I-1)-MUF(I-D) 

TP2  =  DTAN(TETAF(I)+MUF(D) 

F1=-TP2/(TP1-TP2) 

IFCF1.LE.0.)  PRINT  555, I, XI, X2, DX,TP1,TP2, Fl 
555   FORMATC/1X, • I,X1,X2, DX,TP1,TP2, Fl= ' , I5,6D14 .6/) 

IFCF1.LT.0.)  CALL  FINU501) 

DYT=F1*DX*TP1 

IFCDYT.LE.O.)  CALL  FIN(1502) 

DYMIN=DMIN1(DYMIN,STAB*DYT) 
1     CONTINUE 

DYNEXT=DYMIN 

RETURN 

END 


MOVE 


ET1225 
ET1226 
ET1227 
ET1228 
ET1229 
ET1230 
ET1231 
ET1232 
ET1233 
ET1234 
ET1235 
ET1236 
ET1237 
ET1238 
ET1239 
ET1240 
ET1241 
ET1242 
ET1243 
ET1244 
ET1245 
ET1246 
ET1247 
ET1248 
ET1249 
ET1250 
ET1251 
ET1252 
ET1253 
ET1254 
ET1255 
ET1256 
ET1257 
ET1258 
ET1259 
ET1260 
ET1261 
ET1262 
ET1263 
ET1264 
ET1265 


SUBROUTINE  MOVE 

SUBROUTINE  NUMBER   16 

IMPLICIT  REAL*S(A-H,L-Z,$) 

REAL*4  XPL,YPL 

COMMON  /PLUME/XPL(1002,10),YPL(1002) 

COMMON  /IPLUME/KPL,ITYPL(10) 

COMMON  /VECS/XF(101),RMF(101),RPF(101),MF(101),MUF(101), 

1  TETAF(101),BF(101), 

2  XN(101),RMN(101),RPN(101),MN(101),MUN(101), 

3  TETAN(101),BN(101),XTEMP(101) 
COMMON/THICKY/XTH(10  02),TH(1002) 


JET1266 
JET1267 
JET1268 
JET1269 
JET1270 
JET1271 
JET1272 
JET1273 
JET1274 
JET1275 
JET1276 

COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15, JET1277 
G16,G17,G18,G19,G20  JET1278 

JET1279 
JET1280 
JET1281 
JET1282 
JET1283 
JET1284 
JET1285 
JET1286 
JET1287 
JET1288 
JET1289 
JET1290 
JET1291 
JET1292 
JET1293 
JET1294 
JET1295 
JET1296 


G16,G17,G18,G19,G20 
COMMON  /PAR/PAI,PAI2,DEG,XC,YC,MEXIT,MFIN,YMAX,DY0,DY, DYNEXT, 

1  STAB, DEL TA,PSI1,PSIF,ZETA1, SIGMA, FRACG, EPSIL , NUO, 

2  TETSYM,TETLIM,DDY,DYMAX 
/STAG/RHOO,NO,PO,TO,AO,MDOT1 
/IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP,J, 

KF2,IDEL,JDEL,JYXI,JXI,ILEAD,ILEADF,KCLEAD 
/ROW/YF,YN,DXF,DXN 
/CHARAC/XCHARF(92),YCHARF(92),XCHARN(92),YCHARN(92), 

1  RMCARF(92),RPCARF(92),RMCARN(92),RPCARN(92), 

2  TCHARF(92),TCHARN(92),MUCARF(92),MUCARN(92), 

3  CSIGNN(92),CSIGNF(92),MCHARN(92),MCHARF(92), 

4  MCHARK92) 
/ICHARA/KCHARP, KCHARM, KCHARO 
LINE  (N)  IN  OLD  LINE  (F). 


COMMON 
COMMON 

COMMON 
COMMON 


COMMON 
STORE  NEW 
KF  =  KN 
KF2=2*KF 
YF=YN 
DO  1  1=1, KN 
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XF(I)=XN(I)  JET1297 

RMF(I)=RMN(I)  JET1298 

RPF(I)=RPN(I)  JET1299 

MF(I)=MN(I)  JET1300 

MUF(I)=MUN(I)  JET1301 

TETAF(I)=TETAN(I)  JET1302 

BF(I)=BN(I)  JET1303 

1  CONTINUE  JET1304 
DO  2  KC=1,KCHAR0  JET1305 
IFCCSIGNNCKO.EQ.O.)  GO  TO  2  JET1306 
XCHARF(KC)=XCHARN(KC)  JET1307 
YCHARF(KC)=YCHARN(KC)  JET1308 
RMCARF(KC)=RMCARN(KC)  JET1309 
RPCARF(KC)=RPCARN(KC)  JET1310 
TCHARFCKC)=TCHARN(KC)  JET1311 
MUCARF(KC)=MUCARN(KC)  JET1312 
MCHARF(KC)=MCHARN(KC)  JET1313 
CSIGNF(KC)=CSIGNN(KC)  JET1314 

2  CONTINUE  JET1315 
RETURN                                           r\1?AC  X  JET1316 

END Urn-^/N JET1317 

SUBKUUIINb  UPACX JET1318 

C   SUBROUTINE  NUMBER   17  JET1319 

IMPLICIT  REAL*8(A-H,L-Z,$)  JET1320 

REAL**  XPL,YPL  JET1321 

COMMON  /PLUME/XPL(1002,10),YPL(1002)  JET1322 

COMMON  /IPLUME/KPL,ITYPL(10)  JET1323 

COMMON  /VECS/XF(101),RMF(101),RPF(101),MF(101),MUFC101),  JET1324 

1  TETAFC101),BF(101),  JET1325 

2  XN(101),RMN(101),RPN(101),MN(101),MUN(101),  JET1326 

3  TETAN(101),BN(101),XTEMP(101)  JET1327 
COMMON/THICKY/XTHC1 002), TH( 1002)  JET1328 
REAL**  YXI,XI,XIPM,XIGRP,XIAPP,XIF  JET1329 
COMMON  /THICKX/YXI(20),XIC101,20),XIPMC101,20),XIGRP(101,20)  JET1330 

1  ,XIAPPC101,20),XIF(101,20)  JET1331 
COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15,JET1332 

1              G16,G17,G18,G19,G20  JET1333 

COMMON  /PAR/PAI,PAI2,DEG,XC,YC,MEXIT,MFIN,YMAX,DY0,DY,DYNEXT,  JET1334 

1  STAB, DELTA, PSI1 , PSIF, ZETA1 , SIGMA, FRACG, EPSIL , NUO ,  JET1335 

2  TETSYM,TETLIM,DDY,DYMAX  JET1336 
COMMON  /STAG/RHO0,N0,P0,T0,A0,MDOTl  JET1337 
COMMON  /CHARAC/XCHARF(92),YCHARF(92),XCHARN(92),YCHARN(92),  JET1338 

1  RMCARF(92),RPCARF(92),RMCARN(92),RPCARN(92),  JET1339 

2  TCHARF(92),TCHARN(92),MUCARF(92),MUCARN(92),  JET1340 

3  CSIGNN(92),CSIGNF(92),MCHARN(92),MCHARF(92),  JET1341 

4  MCHARK92)  JET1342 
COMMON  /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP,J,  JET1343 

1              KF2,IDEL,JDEL,JYXI,JXI,ILEAD,ILEADF,KCLEAD  JET1344 

COMMON  /ROW/YF,YN,DXF,DXN  JET1345 

C   COMPUTE  X-OPACITY.  JET1346 

C   BEGIN  FROM  LIMITING  CHARACTERISTIC  OF  AN  ASSUMED  P.M.  FAN.  JET1347 

C   XIO   ~   THE  THICKNESS  BETWEEN  THE  LIMITING  CHARACTERISTIC  AND  THE    JET1348 

C            BOUNDARY  CHARACTERISTIC  OF  THE  NUMERICAL  COMPUTATION.  JET1349 

DO  12  1=1, KFO  JET1350 

XIF(I,JXI)=XF(I)  JET1351 

XI(I,JXI)=0.  JET1352 

XIPM(I,JXI)=0.  JET1353 

XIGRP(I,JXI)=0.  JET1354 

XIAPP(I,JXI)=0.  JET1355 

12    CONTINUE  JET1356 

IF(J.EQ.l)  GO  TO  1000  JET1357 

PSILIM=TETLIM  JET1358 

XLIM=XC+(YF-YC)/DTAN(PSILIM)  JET1359 

XBOUND=XF(KF)  JET1360 

KPM=10  JET1361 

DX=(XLIM-XBOUND)/DFLOAT(KPM)  JET1362 

SUM=0.  JET1363 

DO  1  1=1, KPM                                         ,  JET1364 

Xl=XBOUND+DFLOAT(I-l)*DX  JET1365 

X2=X1+DX  JET1366 

PS1=PAI2-DATAN((X1-XC)/(YF-YC))  JET1367 

PS2  =  PAI2-DATAN(CX2-XC)/(YF-YO)  JET1368 
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DO)) 
DO)) 


DO/(MEXITxETA)+DSIN(TETA)/M)/DSIN(PSI) 


Q1=(PS1-PSILIM)/G5 

Q2=(PS2-PSILIM)/G5 

IF(I.EQ.KPM)  Q2=1.D-10 

IFCQ2.LT.O.)  CALL  FINC1701) 

F1=G11X(DSIN(Q1))XX(2.D0/(G-1 

F2=G11X(DSIN(Q2))XX(2.D0/(G-1 

SUM=SUM+DXX(F1+F2)/2.D0 

1  CONTINUE 
XIO=SUMX(NOXSIGMA) 

:   RE-EVALUATE  XIO  FOR  A  RING-JET. 
IFCDELTA.EQ.O.)  GO  TO  14 
M=MFIN 

CALL  MFUNC(M,F,ETA,TETA) 
PSI=TETA+DARSIN(l.DO/M) 
GOREM=1.DO+G1XMXX2 
GOR=MXX2-l.D0 
CALL  HINTER(M,HM) 
DELTOB=0.5DOXDSQRT(GOR)X<1 
1    +G15XHM/2.DO 

EVER=SIGMA*N0*YC/(M*DSINCTETA)*DSIN(PSI)XGOREMXXG6) 
GGG=2.D0-DELT0BX(G+1.D0)/2.D0 

IFCDABS(GGG) .GT.l.D-10)  GO  TO  15 

PRINT  16,  DELTOB,G,GGG 
16    FORMATC/1X, 'FROM  OPACX.  GGG  NEARLY  ZERO.  EXPRESSION  FOR  XIO  IS' 
1         1X,'SINGULAR.   DELT0B,G,GGG=',3D12.4/) 

CALL  FIN(1715) 
15    CONTINUE 

EVER=EVER/GGG 

XI0=EVERX((YF/YC)XXGGG-1.D0)/(YF/YC) 
14    CONTINUE 

XI(KF,JXI)=XIO 

XIPM(KF,JXI)=XIO 

XIGRP(KF,JXI)=XIO 

KF1=KF-1 

DO  2  11=1, KF1 

I=KF-II+1 

X1=XF(I) 

X2=XF(I-1) 

DX=X1-X2 

F1  =  1.D0/(1.D0+G1XMF(I   )XX2)X*G6 

F2=1.D0/(1.D0+G1*MF(I-1)XX2)**G6 
DTNUM=(NOXSIGMA)XDXX(Fl+F2)/2.DO 

XI(I-1,JXI)=XI(I,JXI)+DTNUM 
XIPM(I-1,JXI)=1.D24 
XIGRP(I-1,JXI)=1.D24 
PS1=PAI2-DATAN((X1-XC)/(YF-YC)) 
PS2  =  PAI2-DATAN((X2-XC)/(YF-YO) 
IFCPS2.GT.PSI1)  GO  TO  2 
Q1=(PS1-PSILIM)/G5 
Q2=(PS2-PSILIM)/G5 
IFCQ1.LT.O.)  CALL  FIN(1711) 
F1=G11*(DSIN(Q1))**(2.D0/(G-1.D0)) 
F2=G11*(DSIN(Q2))**(2.D0/(G-1.D0)) 
DTPM=(N0*SIGMA)*DX*(F1+F2)/2.D0 
XIPM(I-1,JXI)=XIPM(I,JXI)+DTPM 
DIST1=DSQRTCCX1-XC)XX2+(YF-YC)**2) 
DIST2=DSQRT((X2-XC)**2+(YF-YC)**2) 
KC1=KCLEAD+I-ILEAD 
KC2=KC1— 1 

IF(KC2.LT.KCLEAD)  GO  TO  21 
M1=MCHARI(KC1) 
M2=MCHARI(KC2) 

CALL  MATCHCI   ,M1 ,MG1 ,MOBIl ,MABI1 ) 
CALL  MATCH(I-l,M2,MG2,MOBI2,MABI2) 
Fl=l.DO/(l.DO+Gl*MGl**2)xxG6 
F2=1.D0/(1.D0+G1*MG2**2)**G6 
DTGRP=(N0*SIGMA)*DX*(F1+F2)/2.D0 
XIGRP(I-1,JXI)=XIGRP(I,JXI)+DTGRP 
21    CONTINUE 

2  CONTINUE 
J   APPROXIMATE  THICKNESS  XIAPPC I , JXI ) .  BASED 

DO  3  1=1, KF 


ON  CLOSED-FORM  INTEGRATION 


JET1369 
JET1370 
JET1371 
JET1372 
JET1373 
JET1374 
JET1375 
JET1376 
JET1377 
JET1378 
JET1379 
JET1380 
JET1381 
JET1382 
JET1383 
JET1384 
JET1385 
JET1386 
JET1387 
JET1388 
JET1389 
JET1390 
JET1391 
JET1392 
JET1393 
JET1394 
JET1395 
JET1396 
JET1397 
JET1398 
JET1399 
JET1400 
JET1401 
JET1402 
JET1403 
JET1404 
JET1405 
JET1406 
JET1407 
JET1408 
JET1409 
JET1410 
JET1411 
JET1412 
JET1413 
JET1414 
JET1415 
JET1416 
JET1417 
JET1418 
JET1419 
JET1420 
JET1421 
JET1422 
JET1423 
JET1424 
JET1425 
JET1426 
JET1427 
JET1428 
JET1429 
JET1430 
JET1431 
JET1432 
JET1433 
JET1434 
JET1435 
JET1436 
JET1437 
JET1433 
.JET1439 
JET1440 
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XIAPPCI,JXI)=1.D24  JET1441 

KC=KCLEAD+(I-ILEAD)  JET1442 

IFCDELTA.EQ.O.)  GO  TO  3  JET1443 

IF(KC.LT.KCLEAD)  GO  TO  3  JET1444 

IF(XFCI).LT.XCHARFCl))  GO  TO  3  JET1445 

M=MCHARI(KC)  JET1446 

CALL  MFUNC(M,F,ETA,TETA)  JET1447 

PSI=TETA+DARSIN(l.DO/M)  JET1448 

GOREM=1.DO+G1*M**2  JET1449 

G0R=M**2-1.D0  JETI450 

CALL  HINTER(M,HM)  JET1451 
DELTOB=0.5DO*DSQRT(GOR)x(l.DO/(MEXIT*ETA)+DSIN(TETA)/M)/DSIN(PSI)  JET1452 

1    +G15XHM/2.DO  JET1453 

EVER=SIGMA*N0*YC/(M*DSIN(TETA)*DSIN(PSI)*GOREM**G6)  JET1454 

GGG=2.D0-DELT0B*(G+1.D0)/2.D0  JET1455 

IFCDABS(GGG).GT.l.D-lO)  GO  TO  25  JET1456 

PRINT  26,  I,KC,M,DELTOB,G,GGG  JET1457 
26    FORMAT(/lX,»FROM  OPACX.  GGG  NEARLY  ZERO.  EXPRESSION  FOR  XIO  IS',   JET1458 

1  IX, 'SINGULAR.   I,KC,M=M5,D12.4/  JET1459 

2  1X,'DELT0B,G,GGG=,,3D12.4/)  JET1460 
CALL  FINC1725)  JET1461 

25    CONTINUE  JET1462 

EVER=EVER/GGG  JET1463 

XIAPP(I,JXI)=EVERX((YF/YC)*XGGG-l.DO)/(YF/YC)  JET1464 

3     CONTINUE  JET1465 

1000  CONTINUE  JET1466 

RETURN                                       /  n/LTsC  JET1467 

END L  UtW*- JET1468 

SUBROUTINE  LOADC  JET1469 

C   SUBROUTINE  NUMBER   18  JET1470 

IMPLICIT  REAL*8CA-H,L-Z,$)  JET1471 

REAL*4  XPL,YPL  JET1472 

COMMON  /PLUME/XPL(1002,10),YPL(1002)  JET1473 

COMMON  /IPLUME/KPL,ITYPLC10)  JET1474 

COMMON  /VECS/XF(101),RMF(101),RPF(101),MF(101),MUF(101),  JET1475 

1  TETAF(101),BF(101),  JET1476 

2  XN(101),RMN(101),RPN(101),MN(101),MUN(101),  JET1477 

3  TETAN(101),BN(101),XTEMP(101)  JET1478 
COMMON/THICKY/XTH(10  02),TH(1002)  JET1479 
REAL*4  YXI,XI,XIPM,XIGRP,XIAPP,XIF  JET1480 
COMMON  /THICKX/YXI(20),XI(101,20),XIPM(101,20),XIGRPC101,20)  JET1481 

1  ,XIAPP(101,20),XIF(101,20)  JET1482 
COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15,JET1483 

1              G16,G17,G18,G19,G20  JET1484 

COMMON  /PAR/PAI,PAI2,DEG,XC,YC,MEXIT,MFIN,YMAX,DY0,DY,DYNEXT,  JET1485 

1  STAB, DELTA, PSI1 , PSIF, ZETA1 , SIGMA, FRACG, EPSIL , NUO ,  JET1486 

2  TETSYM,TETLIM,DDY,DYMAX  JET1487 
COMMON  /STAG/RHO0,N0,P0,T0,A0,MDOTl  JET1488 
COMMON  /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP, J,  JET1489 

1              KF2,IDEL,JDEL,JYXI,JXI,ILEAD,ILEADF,KCLEAD  JET1490 

COMMON  /ROH/YF,YN,DXF,DXN  JET1491 

COMMON  /CHARAC/XCHARF(92),YCHARF(92),XCHARN(92),YCHARN(92),  JET1492 

1  RMCARF(92),RPCARF(92),RMCARN(92),RPCARN(92),  JET1493 

2  TCHARF(92),TCHARN(92),MUCARF(92),MUCARN(92),  JET1494 

3  CSIGNN(92),CSIGNF(92),MCHARN(92),MCHARF(92),  JET1495 

4  MCHARK92)  JET1496 
COMMON  /ICHARA/KCHARP,KCHARM,KCHARO  JET1497 

C   LOAD  FLOW  VARIABLES  OF  GRID  POINTS  IN  THE  FAN  SEGMENT  FROM  THE  JET1498 

C   SEMI-INVERSE  INTEGRATION  (IN  SUBR.  SEMINV).   NOTE  THAT  GRID  POINTS    JET1499 

C   XN(I)  WERE  ALREADY  DETERMINED  IN  SUBR.  GRIDN.  JET1500 

DO  1  I=ILEAD,KN  JET1501 

KC=KCLEAD+I-ILEAD  JET1502 

IF(KC.GT.KCHARP)  CALL  FINU801)  JET1503 

RMN(I)=RMCARN(KC)  JET1504 

RPN(I)=RPCARN(KC)  JET1505 

MN(I)=MCHARN(KC)  JET1506 

MUN(I)=MUCARN(KC)  JET1507 

TETAN(I)=TCHARN(KC)  JET1508 

1     CONTINUE  JET1509 

RETURN                     -                  klUFUtiC  JET1510 

END     N  u  r  u  'v  ^ JET1511 

DOUBLE  PRECISION  FUNCTION  NUFUNCTMl  TETT5T2" 
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SUBROUTINE  NUMBER   19  JET1513 

IMPLICIT  REAL*8CA-H,L-Z,$)  JET1514 

REAL*4  XPL,YPL  JET1515 

COMMON  /PLUME/XPL(1002,10),YPLC1002)  JET1516 

COMMON  /IPLUME/KPL,ITYPL(10)  JET1517 

COMMON  /VECS/XF(101),RMF(101),RPF(101),MF(101),MUF(101),  JET1518 

1  TETAFU01),BFC101),  JET1519 

2  XN(101),RMNU01),RPN(101),MN(101),MUNC101),  JET1520 

3  TETAN(101),BN(101),XTEMP(101)  JET1521 
COMMON/THICKY/XTHC1 002), TH( 1002)  JET1522 
REAL*4  YXI,XI,XIPM,XIGRP,XIAPP,XIF  JET1523 
COMMON  /THICKX/YXI(20),XI(101,20),XIPM(101,20),XIGRP(101,20)  JET1524 

1  ,XIAPP(101,20),XIF(101,20)  JET1525 
COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15,JET1526 

1              G16,G17,G18,G19,G20  JET1527 

COMMON  /PAR/PAI,PAI2,DEG,XC,YC,MEXIT,MFIN,YMAX,DY0,DY,DYNEXT,  JET1528 

1  STAB, DELTA, PSI1 , PSIF,ZETA1 , SIGMA, FRACG, EPSIL , NUO,  JET1529 

2  TETSYM,TETLIM,DDY,DYMAX  JET1530 
COMMON  /STAG/RHO0,N0,P0,T0,A0,MDOTl  JET1531 
COMMON  /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP,J,  JET1532 

1              KF2,IDEL,JDEL,JYXI,JXI,ILEAD,ILEADF,KCLEAD  JET1533 

COMMON  /ROW/YF,YN,DXF,DXN  JET1534 

COMMON  /CHARAC/XCHARF(92),YCHARF(92),XCHARN(92),YCHARN(92),  JET1535 

1  RMCARF(92),RPCARF(92),RMCARN(92),RPCARN(92),  JET1536 

2  TCHARF(92),TCHARN(92),MUCARF(92),MUCARN(92),  JET1537 

3  CSIGNN(92),CSIGNF(92),MCHARN(92),MCHARF(92),  JET1538 

4  MCHARK92)  JET1539 
COMMON  /ICHARA/KCHARP,KCHARM,KCHARO  JET1540 

COMPUTE  NU  AS  FUNCTION  OF  MACH  NUMBER  M.   NOTE  THAT  THE  P.M.  JET1541 
DEFINITION  OF  NU  HAS  BEEN  MODIFIED  BY  ADDING  A  CONSTANT.   THE  USUAL   JET1542 

CHOICE  OF  THE  CONSTANT  IS  SUCH  THAT  NU=0  FOR  INFINITE  M.  JET1543 

Q=1.D0/DSQRT(M**2-1.D0)  JET1544 

NUFUNC=NU0-CG5*DATAN(G5*Q)-DATAN(Q))  JET1545 

RETURN                                          \lkA<ZtrT  JET1546 

END                                               rfNOfc  '  JET1547 


SUBROUTINE  HM5E1 

C   SUBROUTINE  NUMBER   20 

IMPLICIT  REAL*8(A-H,L-Z,$) 
REAL*8  KAPAOB 

COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15, 
1  G16,G17,G18,G19,G20 

COMMON  /PAR/PAI,PAI2,DEG,XC,YC,MEXIT,MFIN,YMAX,DY0,DY,DYNEXT, 

1  STAB, DELTA, PSI1 , PSIF, ZETA1 , SIGMA, FRACG, EPSIL , NUO, 

2  TETSYM,TETLIM,DDY,DYMAX 
COMMON  /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP,J, 

1  KF2,IDEL,JDEL,JYXI,JXI,ILEAD,ILEADF,KCLEAD 

COMMON  /GRP/DMINV,MHINV(101),HMV(101) 
COMMON  /IGRP/KHM 
C   A  ROUTINE  FOR  THE  C+  DERIVATIVE  DUE  TO  RING  SYMMETRY  (GRP) . 
KHM=51 

IFCKHM.GT.101)  CALL  FINC2001) 
MINV0=1.D0/MEXIT 
DMINV=MINVO/DFLOAT(KHM-1) 
M=MEXIT 
SUM=0. 
KHM1=KHM-1 
DO  1  I=1,KHM1 
MF  =  M 

MHINV(I)=MINVO-DFLOATCI-1)*DMINV 
M=1.D0/MHINV(I) 
DM=M~MF 
M1=M-DM 
M2=M-DM/2.D0 
M3  =  M 

CALL  MFUNC(M1,F1,ETA1,TETA1) 
CALL  MFUNC(M2,F2,ETA2,TETA2) 
CALL  MFUNC(M3,F3,ETA3,TETA3) 
SUM=SUM+DM*(Fl+4.D0*F2+F3)/6 .DO 
ETA=ETA3 
TETA=TETA3 
PSI=TETA+DARSIN(1 .DO/M) 

NORM=((3.D0-G)/4.D0)X(M**2-l.D0)**0.75DO/ 


JET1548 
JET1549 
JET1550 
JET1551 
JET1552 
JET1553 
JET1554 
JET1555 
JET1556 
JET1557 
JET1558 
JET1559 
JET1560 
JET1561 
JET1562 
JET1563 
JET1564 
JET1565 
JET1566 
JET1567 
JET1568 
JET1569 
JET1570 
JET1571 
JET1572 
JET1573 
JET1574 
JET1575 
JET1576 
JET1577 
JET1578 
JET1579 
JET1580 
JET1581 
JET1582 
JET1583 
JET1584 
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1     (DSIN(PSI)X(1.D0+G1*MX*2)XXG1«)  JET1585 

HM=SUMXNORM  JET1586 

HMV(I)=HM  JET1587 

GOREM=1.DO+G1XM*X2  JET1588 

GOR=MXX2-1.DO  JET1589 
DELTOB=0.5DOXDSQRT(GOR)X(1.DO/(MEXITXETA)+DSIN(TETA)/M)/DSIN(PSI)  JET1590 

1    +((G+1.D0)/C2.D0X(3.D0-G)))XHM  JET1591 

EPSIOB=DELTOB/DSQRT(GOR)-DSIN(TETA)/(MXDSIN(PSD)  JET1592 

KAPAOB=l.DO  JET1593 

IF(DABS(PAI2-TETA).GT.l.D-6)  JET1594 

1KAPA0B=DTAN(TETA)XEPSI0B  JET1595 

LAMDOB=EPSIOB-DELTOBXGOREM/(GORXDSQRT(GOR))  JET1596 

PRINT  11,I,M,HM,TETAXDEG,PSIXDEG  JET1597 

11  FORMATC/1X,'  I,M,HM,TETA,PSI=',I5,5D12.4)  JET1598 
PRINT  12,DELT0B,EPSI0BXDEG,KAPA0BXDEG,LAMD0BXDEG  JET1599 

12  FORMATC  IX, ' DELTOB, EPSIOB, KAPAOB, LAMD0B= ' , 5X, 5D12 . 4)  JET1600 

I  CONTINUE  JET1601 
MHINV(KHM)=0.  JET1602 
HMV(KHM)=l.DO  JET1603 
RETURN                                       MCI/ Mr  JET1604 

END         Hr  UA/O JET1605 

~~ SUBROUTINE  MFUNC(M,F,  ETA,TETA) JEI16U6 

C   SUBROUTINE  NUMBER   21  JET1607 

IMPLICIT  REALX8(A-H,L-Z,$)  JET1608 

COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15,JET16  09 

1             G16,G17,G18,G19,G20  JET1610 

COMMON  /PAR/PAI,PAI2,DEG,XC,YC,MEXIT,MFIN,YMAX,DY0,DY,DYNEXT,  JET1611 

1  STAB, DELTA, PSI1,PSIF,ZETA1, SIGMA, FRACG, EPSIL , NUO,  JET1612 

2  NUPT1,TETLIM  JET1613 
COMMON  /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP,J,  JET1614 

1             KF2,IDEL,JDEL,JYXI,JXI,ILEAD,ILEADF,KCLEAD  JET1615 

COMMON  /GRP/DMINV,MHINV(101),HMV(101)  JET1616 

C  JET1617 

NU=NUFUNC(M)  JET1618 

TETA=NUFUNC(MEXIT)+PAI2-NU  JET1619 

G0REM=1.D0+G1XMXX2  JET1620 

GOR=MX*2-l.D0  JET1621 

F=(MXX2)X(G0REMXXG13)XDSIN(TETA)/G0RXX1 .25D0  JET1622 

G0REM1=1.D0+G1*MEXITXX2  JET1623 

G0R1=MEXITXX2-1.D0  JET1624 

ETA=((GOREM/GOREM1)XXG14)X((GOR1/GOR)XX0.25DO)  JET1625 

RETURN  JET1626 

END  JET1627 

SUBROUTINE  HINTER(M,H)  JET1628 

C   SUBROUTINE  NUMBER   22  JET1629 

IMPLICIT  REAL*8(A-H,L-Z,$)  JET1630 

COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15,JET1631 

1             G16,G17,G18,G19,G20  JET1632 

COMMON  /PAR/PAI,PAI2,DEG,XC,YC,MEXIT,MFIN,YMAX,DY0,DY,DYNEXT,  J ET16 33 

1  STAB, DELTA, PSI1 , PSIF,ZETA1 , SIGMA, FRACG, EPSIL ,  NUO,  JET1634 

2  TETSYM,TETLIM,DDY,DYMAX  JET1635 
COMMON  /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP,J,  JET1636 

1             KF2,IDEL,JDEL,JYXI,JXI,ILEAD,ILEADF,KCLEAD  JET1637 

COMMON  /GRP/DMINV,MHINV(101),HMV(101)  JET1638 

COMMON  /IGRP/KHM  JET1639 

C  COMPUTE  HCM)  BY  INTERPOLATION  JET1640 

MINV=1.D0/M  JET1641 

I=KHM-IDINT(MINV/DMINV-l.D-9)-l  JET1642 

IF(I.GE.l.AND.I.LT.KHM)  GO  TO  1  JET1643 

PRINT  11,I,KHM,M,MEXIT  JET1644 

II  FORMATC/1X, ,I,KHM,M,MEXIT=,,2I5,2D14.6/)  JET1645 
CALL  FINC2201)  JET1646 

1     CONTINUE  JET1647 

F1=(MINV-MHINV(I+1))/DMINV  JET1648 

F2=1.D0-F1  JET1649 

IFCF1.LT.-1.D-9)  CALL  FIN(2210)  JET1650 

IF(F2.LT.-l.D-9)  CALL  FINC2211)  JET1651 

H=F1*HMV(I)+F2XHMV(I+1)  JET1652 

RETURN                                        LA/L-t—r  l±  JET1653 

END         nfrj    <~+t  JET1654 

SUBROUTINE  MATCH(I,MOB,MAB,M0BI,MABI) JblUbb 

C   SUBROUTINE  NUMBER   23  JET1656 
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IMPLICIT  REAL*8(A-H,L-Z,$)  JET1657 

COMMON  /VECS/XF(101),RMFU01),RPF(101),MF(101),MUFC101),  JET1658 

1  TETAF(101),BF(101),  JET1659 

2  XN(101),RMN(101),RPN(101),MN(101),MUNC101),  JET1660 

3  TETAN(101),BN(101),XTEMPC101)  JET1661 
COMMON  /ROW/YF,YN,DXF,DXN  JET1662 
COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15,JET166  3 

1  G16,G17,G18,G19,G20  JET1664 

COMMON  /PAR/PAI,PAI2,DEG,XC,YC,MEXIT,MFIN,YMAX,DY0,DY,DYNEXT,  JET166  5 

1  STAB, DELTA, PSI1 , PSIF, ZETA1 , SIGMA, FRACG, EPSIL , NUO,  JET1666 

2  TETSYM,TETLIM,DDY,DYMAX  JET1667 
COMMON  /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP,J,  JET1668 

1  KF2,IDEL,JDEL,JYXI,JXI,ILEAD,ILEADF,KCLEAD  JET1669 

COMMON  /GRP/DMINV,MHINV(101),HMV(101)  JET1670 

COMMON  /IGRP/KHM  JET1671 

C  COMPUTE  H(M)  AND  THE  ALFA-DERIVATIVES  JET1672 

M=MOB  JET1673 

CALL  MFUNC(M,F,ETA,TETA)  JET1674 

PSI=TETA+DARSIN(1.D0/M)  JET1675 

CALL  HINTER(M,HM)  JET1676 

GOREM=1.DO+G1XMXX2  JET1677 

GOR=M*X2-1.DO  JET1678 
DELT0B=0.5D0XDSQRT(GOR)x(l.D0/(MEXITxETA)+DSINCTETA)/M)/DSIN(PSI)  JET1679 

1    +G15*HM/2.DO  JET1680 

FOB=(G7XGOREM)XXG2/M  JET1681 

FAB=FOBX(YF/YC)xxDELTOB  JET1682 

CALL  AREAF(FAB,MAB)  JET1683 

C   COMPUTE  MABI  FROM  THE  INVERSE  PROBLEM  SOLUTION  JET1684 

COTAV=(XF(I)-XC)/(YF-YC)  JET1685 

PSI0=PAI2-DATAN(C0TAV)  JET1686 

EVY=YFXDLOG(YF/YC)/(YF-YC)-l.D0  JET1687 

PSIN=PSIO  JET1688 

DO  1  ITER=1,50  JET1689 

PSI=PSIN  JET1690 

M=DSQRT(1.D0+G4/DTANUPSI-TETLIM)/G5)XX2)  JET1691 

M=DMAX1(M,MEXIT)  JET1692 

CALL  HINTER(M,HM)  JET1693 

CALL  MFUNC(M,F,ETA,TETA)  JET1694 

GOREM=1.DO+G1XMXX2  JET1695 

G0R=MXX2-1.D0  JET1696 
DELTOB=0.5DOXDSQRT(GOR)*(1.DO/(MEXITXETA)+DSINCTETA)/M)/DSINCPSI)  JET16  97 

1    +G15XHM/2.D0  JET1698 

EPSIOB=DELTOB/DSQRT(GOR)-DSIN(TETA)/(MXDSIN(PSD)  JET16  99 

LAMDOB=EPSI OB-DEL TOBXGOREM/(GORXDSQRT(GOR))  JET17  0  0 

COTN=COTAV+LAMDOBXEVY/DSIN(PSI)XX2  JET17  01 

PSIN=PAI2-DATAN(C0TN)  JET1702 

DPSI=PSIN-PSI  JET1703 

IF(DABS(DPSI).LT.l.D-9)  GO  TO  11  JET1704 

I  CONTINUE  JET1705 
PRINT  12,I,ITER,PSI,PSIN,DPSI,M,XF(I),YF,XC,YC  JET1706 

12    F0RMAT(/1X, • I , ITER, PSI , PSIN, DPSI ,M, XF( I ) , YF, XC, YC= »//  JET1707 

1        1X,2I4,8D11.3/)  JET1708 

CALL  FINC2301)  JET1709 

II  CONTINUE  JET1710 
C   USING  M0BI=M  AS  COMPUTED  FROM  THE  INVERSE  PROBLEM,  FIND  MABI.  JET1711 

M0BI=M  JET1712 

M=M0BI  JET1713 

CALL  MFUNC(M,F,ETA,TETA)  JET1714 

PSI=TETA+DARSIN(1 .DO/M)  JET1715 

CALL  HINTER(M,HM)  JET1716 

GOREM=1.DO+G1XMXX2  JET1717 

G0R=MXX2-1 .DO  JET1713 
DELT0B=0.5D0XDSQRT(GOR)X(l.D0/(MEXITXETA)+DSIN(TETA)/M)/DSIN(PSI)  JET1719 

1    +G15XHM/2.D0  JET1720 

F0B=(G7XGOREM)xxG2/M  JET1721 

FAB=FOBX(YF/YC)xxDELTOB  JET1722 

CALL  AREAF(FAB,MABI)  JET1723 

RETURN  A0P/VP  JET1724 

ENID LL_— JET1725 

SUBROUTINE  AREAF(F,M)  Jtll/26 

C   SUBROUTINE  NUMBER   24  JET1727 

IMPLICIT  REALX8(A-H,L-Z,$)  JET1728 
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COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15,JET1729 

1              G16,G17,G18,G19,G20  JET1730 

COMMON  /PAR/PAI,PAI2,DEG,XC,YC,MEXIT,MFIN,YMAX,DYO,DY,DYNEXT,  JET1731 

1  STAB, DELTA, PSI1 ,PSIF,ZETA1, SIGMA, FRACG, EPSIL, NUO,  JET1732 

2  TETSYM,TETLIM,DDY,DYMAX  JET1733 
COMMON  /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP,J,  JET1734 

1             KF2,IDEL,JDEL,JYXI,JXI,ILEAD,ILEADF,KCLEAD  JET1735 

COMMON  /GRP/DMINV,MHINV(101),HMV(101)  JET1736 

COMMON  /IGRP/KHM  JET1737 

C   COMPUTE  MACH  NUMBER  M  FROM  AREA  RATIO  FUNCTION  F  JET1738 

C   F=C(2/CG+1))X(1+(G-1)*M*X2))*X((G+1)/(2*(G-1)))/M  JET1739 

C   INITIAL  GUESS  IS  MIN  JET1740 

E1=(F*MEXIT)*X(1.D0/G2)/G7  JET1741 

E2=(E1-1.D0)/G1  JET1742 

E3=DMAX1(E2,MEXIT*X2)  JET1743 

MIN=DSQRT(E3)  JET1744 

EMN=MIN  JET1745 

DO  1  1=1,100  JET1746 

EMO  =  EMN  JET1747 

GOREM=1.DO+G1*EMOX*2  JET1748 

GOR=EMO**2-1.DO  JET1749 

F0=(G7*G0REM)**G2/EM0  JET1750 

DF=FO-F  JET1751 

C     PRINT  123,I,EM0,EMN,F0,F,DF,G0R,G0REM  JET1752 

C123   F0RMATC1X, • I , EMO, EMN, FO, F, DF, GOR, GOREM= ' , I5,7D12 . 4)  JET1753 

DFDM=FOXGOR/(EMO*GOREM)  JET1754 

DMN=DF/DFDM  JET1755 

EMN=EMO-DMN  JET1756 

EPSEM=DABS(DMN/EMN)  JET1757 

IFCEPSEM.LT. l.D-10)  GO  TO  11  JET1758 

I  CONTINUE  JET1759 
CALL  FINC2401)  JET1760 

II  CONTINUE  JET1761 
M=EMN  JET1762 
RETURN  JET1763 
END  JET1764 

■ 
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